<figure><img src="image_1.png" /></figure>

HAL
open science

# Flutter analysis of fixed and rotary wings through a one-dimensional unified formulation

M. Filippi, E. Carrera

► To cite this version:

M. Filippi, E. Carrera. Flutter analysis of fixed and rotary wings through a one-dimensional unified formulation. Composite Structures, 2015, 133, pp.381 - 389. 10.1016/j.compstruct.2015.07.103. hal-01398281

HAL Id: hal-01398281

https://hal.science/hal-01398281

Submitted on 17 Nov 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Flutter analysis of fixed and rotary wings through a one-dimensional
unified formulation

M. Filippi *, E. Carrera

Department of Mechanical and Aerospace Engineering, Politecnico di Torino, Corso Duca degli Abruzzi 24, 10129 Torino, Italy

Carrera Unified Formulation (CUF) is used to perform flutter analyses of fixed and rotary wings. The one dimensional refined theories are obtained through an axiomatic enrichment of the displacement field components by only setting the input parameters, namely the number of terms and the kind of the cross sectional functions. Within this work, Taylor like expansions of N order (TEN) are used. The aerodynamic loadings are determined through the unsteady strip theories proposed by Theodorsen and Loewy. The finite element method is used to solve the governing equations that are derived, in a weak form, using the generalized Hamilton's Principle. These equations are written in terms of CUF "fundamental nuclei", which do not vary with the theory order (N). The flutter instability of fixed and rotary wings with rectangular and realistic cross sections is investigated. The results are reported in terms of flutter velocities and frequencies and, when possible, they are compared with experimental, numerical and analytical solutions. Despite the intrinsic limitations of the used aerodynamic theories, the proposed methodology appears valid for aeroelastic and vibrational analyses of several structures by ensuring a significant accuracy with a low computational cost.

**Keywords:** Rotating beams, Finite element method, Higher-order theories, Composite structures

## 1. Introduction

According to Collar's definition, aeroelasticity is “the study of the mutual interaction that takes place within the triangle of the inertial, elastic, and aerodynamic forces acting on structural members exposed to an air stream, and the influence of this study on design”. The fluid structure interaction (FSI) has brought to catastrophic events due to sudden failures of bridges, airplanes, helicopter blades, etc. A correct and safe design must require, therefore, an accurate prediction of aeroelastic phenomena. Unfortunately, FSI analyses are, in most cases, too computationally expensive; thus a tradeoff between accuracy and cost is often necessary.

For fixed wing configurations, the simplest aerodynamic theories based on strip approaches were proposed in the first half of the twentieth century [1, 5]. These models were properly combined with simplified structural theories to develop reliable aeroelastic tools. For instance, aeroelastic analyses were performed on wings [6, 8] and civil structures [9] by using the beam plate

approach. On the basis of this methodology, the flexural, torsional and secondary stiffness of the cross section are reduced to the equivalent beam quantities (EI, GJ, K, etc.). To overcome the intrinsic limitations of this technique (no coupling due to the material distribution), Librescu et al. proposed a refined 1D theory including non classical effects, such as warping and transverse shear defor mations and the non uniform torsion [10 12]. Moreover, Patil et al. developed a non linear 1D formulation for laminated box beams by adopting an asymptotic technique to compute the cross sectional mechanical properties [13,14]. Later, Palacios et al. used the same formulation to evaluate the aeroelastic behav ior of highly flexible wings [15]. Nowadays, the rapid increase of the computational power has motivated the development in the area of computational fluid dynamics (CFD). Over the last 30 years, the CFD based aeroelasticity progressed from full potential theo ries (strip theory, panel methods) to problems governed by the Navier Stokes equations. Recent works offer interesting overviews of CFD aeroelastic tools [16,17] furthermore providing useful information about the emerging trends in the aeroelasticity field [18]. From a structural point of view, the CFD simulations can be coupled with finite element models, in which non linear elements may be used to include either geometrical effects or large deforma tions. A detailed paper about non linear structural models can be found in [19].
Regarding the aeroelasticity of rotary wings, the problem becomes more complex due to the geometric nonlinearities that must be taken into account in both elastic and aerodynamic terms. This need has driven the development of suitable structural and aerodynamic theories as confirmed by the considerable number of articles available in the literature. Almost 50 years ago, Loewy provided a thorough overview on a broad range of topics related to the dynamics and aeroelasticity of rotary wings such as flap lag flutter, pitch lag flutter and ground and air resonance [20]. Although being more limited in scope, other contributions to the study of classical flutter conditions and unsteady aerodynamic theories were proposed in [21,22], respectively. In that period, the aerodynamic models for the instability predictions were essentially extensions of steady and quasi steady strip theories conceived for fixed wing configurations. The first important unsteady approach for hover, based on Theodorsen's model, was proposed by Loewy [23]. Though in an approximate manner, this theory takes into account the effect of the spiral returning wake beneath the rotor. Other valuable models were proposed on the basis of Greenberg's theory [24], where a pulsating velocity variation and a constant pitch angle were included in Theodorsen's model. These modifications were extended to the case of rotary wings in [25], in which Theodorsen's, Loewy's and Poisso's theories were used for studying the flap lag torsional coupling. The authors pointed out that the above modifications should be included to realistically reproduce the aeroelasticity of rotary wings, where the assumptions commonly used in deriving strip theories (1 cross sections are assumed to perform only simple harmonic pitching and plunging oscillations about a zero equilibrium position; 2 the velocity of oncoming airflow is constant; 3 usual potential small disturbance unsteady aerodynamics are assumed to apply) are intrinsically violated. The modified strip theory with the quasi steady approximation was adopted to study the stability of composite helicopter blades with swept tips both in hover and in forward flight [26]. The swept tip blades were modeled using non linear 1D finite element with the inclusion of transverse shear deformations and out of plane warping. These results were recently used to verify the accuracy of a nonlinear structural for mulation valid for slender, homogeneous and twisted blades [27]. Other analyses on advanced geometry blades were carried out in [28] with the purpose of evaluating the effects of sweep and droop on both rotor aeroelastic stability and rotorcraft aero mechanical stability. The evolution of the mechanics of helicopter blades, focusing on the aeroelastic and aerodynamic issues, was extensively discussed in [29, 32]. Moreover, even if the aeroelasticity of large wind turbine blades is inherently different, several related aeroelastic problems were solved using rotary wing theories [33, 35]. Furthermore, the response problems of an isolated wind turbine blade and a complete rotor/tower configuration were detailed discussed in [36, 38].

Within this work, we propose an advanced 1D formulation to predict the flutter conditions of various fixed and rotary wing con figurations. The higher order beam theories are obtained using the Carrera Unified Formulation (CUF), which enables to derive, at least theoretically, an infinite number of models [39, 41]. The Taylor type expansions (referred to as TE) represent a particular class of the 1D CUF theories. The TE models have been adopted for the study of mechanical behaviors of thin walled structures [42,43] and non conventional cross sections made of isotropic, composite [39,40] and functionally graded [41] materials. Other encouraging results have also been obtained in the dynamics of rotors by analyzing both blades [44] and spinning shafts [45, 47]. Furthermore, aeroelastic analyzes were performed on fixed wing configurations by combining the TE elements with aerodynamic panel methods [48,49] and unsteady strip theories [50].

## 2. The structural model: Carrera Unified Formulation

CUF states that the displacement field, $\mathbf{u}^T(x, y, z, t)$ ($u_x, u_y, u_z)^T$, is an expansion of generic functions, $F_τ(x, z)$ for the displacement vector, $\mathbf{u}_τ(y)$:

$$ \mathbf{u}(x, y, z, t) \quad F_{\tau}(x, z) \mathbf{u}_{\tau}(y, t) \quad \tau \quad 1, 2, \dots, T \quad (1) $$

where T is the number of terms in the expansion and, according to Einstein's generalized notation, $τ$ indicates summation. In this work, Eq. (1) consists of Taylor like expansions denoted as 'TE'. For example, the third order displacement field (TE3) is:

$$ u_x = u_{x1} + x u_{x2} + z u_{x3} + x^2 u_{x4} + xz u_{x5} + z^2 u_{x6} + x^3 u_{x7} + x^2 z u_{x8} + xz^2 u_{x9} + z^3 u_{x10} $$

$$ u_y = u_{y1} + x u_{y2} + z u_{y3} + x^2 u_{y4} + xz u_{y5} + z^2 u_{y6} + x^3 u_{y7} + x^2 z u_{y8} + xz^2 u_{y9} + z^3 u_{y10} $$

$$ u_z = u_{z1} + x u_{z2} + z u_{z3} + x^2 u_{z4} + xz u_{z5} + z^2 u_{z6} + x^3 u_{z7} + x^2 z u_{z8} + xz^2 u_{z9} + z^3 u_{z10} $$

The remarkable feature of these models is that classical beam theories can be obtained as particular cases of Taylor expansions. Although, the capabilities of the TE approaches have been evaluated in several works, we suggest referring to [51] for a detailed description. The stresses and the strains are grouped as follows:

$$ \begin{aligned} \epsilon_p &= \{\epsilon_{zz} \quad \epsilon_{xx} \quad \epsilon_{xz}\}^T \quad \sigma_p = \{\sigma_{zz} \quad \sigma_{xx} \quad \sigma_{xz}\}^T \\ \epsilon_n &= \{\epsilon_{zy} \quad \epsilon_{xy} \quad \epsilon_{yy}\}^T \quad \sigma_n = \{\sigma_{zy} \quad \sigma_{xy} \quad \sigma_{yy}\}^T \end{aligned} \quad (2) $$

where the subscript "p" stands for the terms lying on the cross section, while "n" stands for the terms lying on the other planes, which are orthogonal to the cross section. The strain displacement relations and Hooke's law are, respectively:

$$ \begin{aligned} \epsilon_p &= \mathbf{D}_p \mathbf{u} \\ \epsilon_n &= \mathbf{D}_n \mathbf{u} \quad (\mathbf{D}_{np} + \mathbf{D}_{ny})\mathbf{u} \end{aligned} \quad (3) $$

$$ \begin{aligned} \sigma_p &= \tilde{C}_{pp} \epsilon_p + \tilde{C}_{pn} \epsilon_n \\ \sigma_n &= \tilde{C}_{np} \epsilon_p + \tilde{C}_{nn} \epsilon_n \end{aligned} \quad (4) $$

where $D_p$ and $D_n$ are linear differential operators. Both laminated box beams and cylinders can be considered constituted by a certain number of either straight or curved plates of orthotropic material, whose material coordinate systems generally do not coincide with the physical coordinate system $(x,y,z)$ of the structure. The matrices of material coefficients of the generic material k are

$$ \tilde{C}_{pp}^k = \begin{bmatrix} \tilde{C}_{11}^k & \tilde{C}_{12}^k & \tilde{C}_{14}^k \\ \tilde{C}_{12}^k & \tilde{C}_{22}^k & \tilde{C}_{24}^k \\ \tilde{C}_{14}^k & \tilde{C}_{24}^k & \tilde{C}_{44}^k \end{bmatrix}, \quad \tilde{C}_{pn}^k = \begin{bmatrix} \tilde{C}_{15}^k & \tilde{C}_{16}^k & \tilde{C}_{13}^k \\ \tilde{C}_{25}^k & \tilde{C}_{26}^k & \tilde{C}_{23}^k \\ \tilde{C}_{45}^k & \tilde{C}_{46}^k & \tilde{C}_{43}^k \end{bmatrix}, \quad (5) $$

For the sake of the brevity, the explicit forms of the coefficients of the C matrices are shown in [46]. The finite element method is used in order to consider arbitrary cross section profiles. The generalized displacement vector is

$$ \mathbf{u}_{\tau}(y, t) \quad N_i(y) \mathbf{q}_{\tau i}(t) \quad (6) $$

where $N_i(y)$ are the shape functions and $\mathbf{q}_{\tau i}(t)$ is the nodal displacement vector

$$ \mathbf{q}_{\tau i}(t) = \{ q_{u_{x\tau i}} \quad q_{u_{y\tau i}} \quad q_{u_{z\tau i}} \}^T \quad (7) $$
### 3. 2-D steady and unsteady aerodynamic theories for fixed-wing

Despite steady and unsteady strip theories were already combined with the 1D CUF elements in [50], for the sake of complete ness, this section provides a brief description of these aerodynamic approaches. According to the assumptions of the small disturbance theory, Theodorsen considered a flat plate with a control surface, which moves in vertical translation $h(t)$ and rotate about an axis at $x$ $ba$ through an angle $\alpha(t)$. As derived in [52], the final expression of the lift for unit of span is

$$L \quad \pi \rho_a b^2 [\ddot{h} + U\dot{\alpha} \quad ba\ddot{\alpha}] + 2\pi \rho_a UbC(k) \left[\dot{h} + U\alpha + b\left(\frac{1}{2} \quad a\right)\dot{\alpha}\right] \quad (8)$$

where 'U' is the free stream velocity, 'b' the semihorod, '$\rho_a$' the air density and '$C(k)$' the deficiency lift function, which depends on the so called reduced frequency $k$ $\frac{\omega U}{b}$. The term 'a' defines the position of the rotation axis with respect to the center of the section and it is dependent on the support condition, the lamination scheme and the applied load. Due to the difficulty involved in defining its correct value for swept and asymmetric laminated structures, the elastic offset 'a' is assumed to be null. Theodorsen identified the term '$C(k)$' as a Hankel function of the second kind, which is, in turn, comprised of Bessel functions of the first and second kinds [52,53]

$$C(k) \quad F(k) + iG(k) \quad \frac{H_1^{(2)}(k)}{H_1^{(2)}(k) + iH_0^{(2)}(k)} \quad (9)$$

The first contribution of Eq. (8), contains the virtual mass (or non circulatory terms) that acts as a form of inertia of the fluid surrounding the airfoil and, since the single and doubly differentiated terms relative to the mass properties of the structure are small, they are usually neglected. The remaining terms determine the value of the lift and represent the circulatory part. The steady state simplification can be easily obtained by replacing the function $C(k)$ with the unitary value and omitting the term containing $\dot{\alpha}$.

$$L \approx 2\pi \rho_a Ub[\dot{h} + U\alpha] \quad (10)$$

Diederich's approximation is used to correct the sectional lift coefficient '$C_l$' relating to the influence of the aspect ratio $(AR \frac{2L_\alpha}{c_m})$ and the sweep angle ('A')

$$c_{l\alpha} \quad \frac{dC_L}{d\alpha} \quad \frac{\pi AR}{\pi AR + c_{l\alpha 0} \cos(\Lambda)} \quad c_{l\alpha 0} \cos(\Lambda) \quad (11)$$

where '$L_w$' is the length of the wing, '$c$' is the mean chord and, '$c_{l\alpha 0}$' is the lift curve slope equal to $2\pi$. Finally, since the quantity '$b\pi$' may be approximated with the integral $\int_b^b \frac{b \frac{x}{b+x} dx}{dx}$, which reproduces the distribution of pressure over a slightly inclined, thin and uncambered airfoil, Eqs. (8) and (10) become

$$L \quad \frac{2\pi AR \cos(\Lambda)}{\pi AR + 2\pi \cos(\Lambda)} \int_b^b \times \frac{b}{b+x} \frac{x}{dx} \rho_a UC(k) \left[ \dot{h} + U\alpha + b\left(\frac{1}{2} \quad a\right) \dot{\alpha} \right] \quad (12)$$

$$L \approx \frac{2\pi AR \cos(\Lambda)}{\pi AR + 2\pi \cos(\Lambda)} \int_b^b \frac{b}{b+x} \frac{x}{dx} \rho_a U[\dot{h} + U\alpha] \quad (13)$$

#### 3.1. Unsteady aerodynamic theory for blade

According to [54], the general expression for the lift distribution on a rotary wing yields

$$\frac{dL}{dy} \quad \pi \rho_a b^2 \Omega^2 \left[ \frac{\ddot{h}}{\Omega^2} + \frac{\dot{\alpha}}{\Omega} \quad ba\frac{\ddot{\alpha}}{\Omega} \right] + 2\pi \rho_a b \Omega^2 C(k) \left[ \frac{\dot{h}}{\Omega} + y\alpha + b\left(\frac{1}{2} \quad a\right) \frac{\dot{\alpha}}{\Omega} \right] \quad (14)$$

where the reduced frequency 'k' is $\frac{2\omega\zeta}{3\Omega R}$ and, 'R' is rotor radius. Although this formulation is limited as it relies on Theodorsen problem, it can be easily generalized since the unsteadiness is described by the lift deficiency function.

In the general case, the rotor blade develops oscillating loads, which generate a radial distribution of vorticity blown downward at velocity 'u'. The wake can be described via models with increasing levels of complexity but, historically, the first successful attempt to incorporate the "returning wake" into the unsteady aerodynamics problem was proposed by Loewy [23]. According to Loewy's theory, the returning wake is modeled as an infinite number of layers, which are equally spaced by a fixed distance, '$h_i$' semi chords. Here, the factor '$h_i$' quantifies the contribution of the underlying vortices to the lift distribution on the reference blade. The main difference with respect to Theodorsen's model lies in the form of the lift deficiency function

$$C'(k, m, h_l) \quad \frac{H_1^{(2)}(k) + 2J_1(k) W(km, h_l)}{H_1^{(2)}(k) + iH_0^{(2)}(k) + 2[J_1(k) + iJ_0(k)]W(km, h_l)} \quad (15)$$

where '$J_n(k)$' and '$H_n^{(2)}(k)$' are the first kind Bessel and second kind Hankel functions of order $n$, respectively. The wake weighting function '$W(k,m,h_l)$' for a rotor with '$Q$' blades is defined as

$$W(k, m, h) \quad \frac{1 + \sum_q^Q \frac{1}{q} (e^{kh_l Q} e^{i2\pi m q})^q e^{i\Phi_q}}{e^{kh_l Q} e^{i2\pi m} \quad 1} \quad (16)$$

where '$m$' and '$\Phi_q$' are the ratio of oscillatory frequency to rotational frequency ($\omega/\Omega$) and the phase angle between the $q$th blade and the reference blade. Due to periodicity, Eq. (16) can be expressed following the formulation for a single blade rotor with the modified values '$\dot{m}$' and '$\dot{h}$'

$$W(k, \dot{m}, \dot{h}) \quad \frac{1}{e^{khQ} e^{i2\pi \dot{m}}} \quad \dot{m} \quad \frac{\omega}{Q\Omega} \quad \dot{h} \quad \frac{2\pi \lambda_0}{Q(b/r)} \quad (17)$$

where '$\lambda_0 \approx \sqrt{C_T/2}$' (where $C_T$ is the thrust coefficient) is the inflow ratio and, '$r$' is the generic radial distance.

### 4. Equations of motion

The general form of Hamilton's Principle states that the variation of the kinetic and potential energy and the variation of the work exerted by non conservative forces must equal zero during any time interval. It may be expressed as

$$\int_{t_0}^{t_1} \delta(T (U + U_{\sigma 0})) dt + \int_{t_0}^{t_1} \delta L_{\text{ext}} dt \quad 0 \quad (18)$$

where $\delta$ refers to kinematically admissible perturbations during the indicated time interval and

* T: total kinetic energy of the system

* $U + U_{\sigma 0}$: total potential energy of the system, including both strain energy ('U') and potential of any conservative external forces ('U_{\sigma 0}')

* $L_{\text{ext}}$: work done by non conservative forces acting on the system, including any arbitrary external loads

The application of this variational concept leads directly to the equations of motion (EoM) for any given system. The energy expressions for a rotary wing are
$$ T \frac{1}{2} \int_V \rho (\dot{\mathbf{u}}^T \dot{\mathbf{u}} + 2\mathbf{u}^T \Omega^T \dot{\mathbf{u}} + \mathbf{u}^T \Omega^T \Omega \mathbf{u} + 2\dot{\mathbf{u}}^T \Omega \mathbf{r} + 2\mathbf{u}^T \Omega^T \Omega \mathbf{r}) dV \quad (19) $$

$$ U \frac{1}{2} \int_V (\mathbf{u}^T \mathbf{D}^T \mathbf{C} \mathbf{D} \mathbf{u}) dV \quad (20) $$

$$ U_{\sigma_0} \int_V (\boldsymbol{\sigma}_0^T \boldsymbol{\epsilon}_{nl}) dV \quad \sigma_0 \quad \Omega^2 \rho \left[ r_h R + \frac{1}{2} R^2 \quad r_h r \quad \frac{1}{2} r^2 \right] \quad (21) $$

where 'ρ' is the material density, 'r' is the distance from the neutral axis, 'r_h' is the hub radius and the constant rotational speed matrix is

$$ \Omega \begin{bmatrix} 0 & \Omega & 0 \\ \Omega & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix} $$

The external work generated by the lift can be written as

$$ L_{\text{ext}} \int_{y} \int_{x} u_z(x, y, z_{\text{top}}) L(x, y, z_{\text{top}}) dx dy \quad (22) $$

where 'ztop' is the upper z coordinate of the cross section. According to Eqs. (8) and (14), Eq. (22) for fixed and rotary wings becomes, respectively

$$ L_{\text{ext}} \quad \text{cost} \left[ \int_{y} \int_{x} \mathbf{u}^{T} \sqrt{\frac{b}{b+x}} \mathbf{I}_{L} \dot{\mathbf{u}} + \mathbf{u}^{T} \sqrt{\frac{b}{b+x}} \mathbf{I}_{L} \mathbf{u}_{x} \right. \\ \left. + \mathbf{u}^{T} b \left(\frac{1}{2} \quad a\right) \sqrt{\frac{b}{b+x}} \mathbf{I}_{L} \dot{\mathbf{u}}_{x} \right] dx dy $$

$$ L_{\text{ext}} \quad \text{cost} \left[ \int_{y} \int_{x} \mathbf{u}^{T} \Omega y \sqrt{\frac{b}{b+x}} \mathbf{I}_{L} \dot{\mathbf{u}} + \mathbf{u}^{T} \Omega^{2} y^{2} \sqrt{\frac{b}{b+x}} \mathbf{I}_{L} \mathbf{u}_{x} dx dy \right] \\ + \text{cost} \left[ \int_{y} \int_{x} \mathbf{u}^{T} \Omega y \sqrt{\frac{b}{b+x}} \mathbf{I}_{L} b \left(\frac{1}{2} \quad a\right) \dot{\mathbf{u}}_{x} dx dy \right] $$

where

$$ \dot{h} \quad \dot{u}_z \quad \propto \quad \frac{du_z}{dx} \quad u_{z,x} \quad \propto \quad \dot{u}_{z,x} $$

$$ \mathbf{I}_L \begin{bmatrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 1 \end{bmatrix} \quad \text{cost} \quad \frac{2\pi A R \cos(\Lambda)}{\pi A R + 2\pi \cos(\Lambda)} \rho_a U C(k) $$

Introducing CUF (Eq. (1)) and FE (Eq. (6)) approximations and substituting Eqs. (19) (22) in Eq. (18), we obtain

$$ \int_{t_0}^{t_1} \left( \delta \mathbf{q}_{\tau i}^T \mathbf{M}^{ij\tau s} \dot{\mathbf{q}}_{sj} + \delta \mathbf{q}_{\tau i}^T \mathbf{D}_{\text{TOT}}^{ij\tau s} \dot{\mathbf{q}}_{sj} + \delta \mathbf{q}_{\tau i}^T \mathbf{K}_{\text{TOT}}^{ij\tau s} \mathbf{q}_{sj} + \delta \mathbf{q}_{\tau i}^T \mathbf{F}_{\Omega}^{\tau s} \mathbf{r} \right) dt \quad 0 \quad (23) $$

where the matrices are written in terms of CUF fundamental nuclei. For fixed wings, we have

$$ \mathbf{D}_{\text{TOT}}^{ij\tau s} \quad \mathbf{D}_{L}^{ij\tau s} + \mathbf{D}_{Lm}^{ij\tau s} $$

$$ \mathbf{K}_{\text{TOT}}^{ij\tau s} \quad \mathbf{K}^{ij\tau s} + \mathbf{K}_{L}^{ij\tau s} $$

whereas, for rotary wings

$$ \mathbf{D}_{\text{TOT}}^{ij\tau s} \quad \mathbf{D}_{L\Omega}^{ij\tau s} + \mathbf{D}_{Lm\Omega}^{ij\tau s} + \mathbf{D}_{\Omega}^{ij\tau s} $$

$$ \mathbf{K}_{\text{TOT}}^{ij\tau s} \quad \mathbf{K}^{ij\tau s} + \mathbf{K}_{L\Omega}^{ij\tau s} + \mathbf{K}_{\Omega}^{ij\tau s} + \mathbf{K}_{\sigma_0}^{ij\tau s} $$

The subscripts 'L' and 'Ω' denote the contributions due to the aero dynamic load and rotation, respectively. The compact expressions of the fundamental nuclei yield

$$ \mathbf{M}^{ij\tau s} \quad I_{l}^{ij} \triangleleft (F_{\tau} \rho \mathbf{I}_{F_{s}}) \triangleright $$

$$ \mathbf{D}_{\Omega}^{ij\tau s} \quad I_{l}^{ij} \triangleleft (F_{\tau} \rho \mathbf{I}_{F_{s}}) \triangleright 2\Omega $$

$$ \mathbf{D}_{L}^{ij\tau s} \quad \text{cost} \ I_{l}^{ij} \triangleleft \left( \sqrt{\frac{b}{b+x}} \frac{x}{F_{\tau}(x, z_{\text{top}})} \mathbf{I}_{L} F_{s}(x, z_{\text{top}}) \right) \triangleright x $$

$$ \mathbf{D}_{Lm}^{ij\tau s} \quad \text{cost} \ I_{l}^{ij} b \left(\frac{1}{2} \quad a\right) \triangleleft \left( \sqrt{\frac{b}{b+x}} \frac{x}{F_{\tau,x}(x, z_{\text{top}})} \mathbf{I}_{L} F_{s}(x, z_{\text{top}}) \right) \triangleright x $$

$$ \mathbf{D}_{L\Omega}^{ij\tau s} \quad \text{cost} \ \Omega I_{ly}^{ij} \triangleleft \left( \sqrt{\frac{b}{b+x}} \frac{x}{F_{\tau}(x, z_{\text{top}})} \mathbf{I}_{L} F_{s}(x, z_{\text{top}}) \right) \triangleright x $$

$$ \mathbf{D}_{Lm\Omega}^{ij\tau s} \quad \text{cost} \ \Omega I_{ly}^{ij} b \left(\frac{1}{2} \quad a\right) \triangleleft \left( \sqrt{\frac{b}{b+x}} \frac{x}{F_{\tau,x}(x, z_{\text{top}})} \mathbf{I}_{L} F_{s}(x, z_{\text{top}}) \right) \triangleright x $$

$$ \mathbf{K}_{\Omega}^{ij\tau s} \quad I_{l}^{ij} \triangleleft (F_{\tau} \rho \mathbf{I}_{F_{s}}) \triangleright \Omega^{T} \Omega $$

$$ \mathbf{K}_{\sigma_0}^{ij\tau s} \quad I_{l\sigma_0}^{ij,y} \triangleleft (F_{\tau} \rho \mathbf{I}_{F_{s}}) \triangleright \Omega^{T} \Omega $$

$$ \mathbf{K}^{ij\tau s} \quad I_{l}^{ij} \triangleleft \mathbf{D}_{np}^{T} (F_{\tau} \mathbf{I}) [\tilde{\mathbf{C}}_{np}^{k} \mathbf{D}_{p}(F_{s} \mathbf{I}) + \tilde{\mathbf{C}}_{nn}^{k} \mathbf{D}_{np}(F_{s} \mathbf{I})] \\ + \mathbf{D}_{p}^{T} (F_{\tau} \mathbf{I}) [\tilde{\mathbf{C}}_{pp}^{k} \mathbf{D}_{p}(F_{s} \mathbf{I}) + \tilde{\mathbf{C}}_{pn}^{k} \mathbf{D}_{np}(F_{s} \mathbf{I})] \triangleright \\ + I_{l}^{ij,y} \triangleleft [\mathbf{D}_{np}^{T} (F_{\tau} \mathbf{I}) + \mathbf{D}_{p}^{T} (F_{\tau} \mathbf{I}) \tilde{\mathbf{C}}_{pn}^{k}] F_{s} \triangleright + \mathbf{I}_{Ay} \\ + I_{l}^{ij,y} \tilde{\mathbf{I}}_{Ay}^{T} \triangleleft F_{\tau} \tilde{\mathbf{C}}_{np}^{k} \mathbf{D}_{p}(F_{s} \mathbf{I}) + \tilde{\mathbf{C}}_{nn}^{k} \mathbf{D}_{np}(F_{s} \mathbf{I}) \triangleright \\ + I_{l}^{ij,y} \tilde{\mathbf{I}}_{Ay}^{T} \mathbf{I}_{Ay} \triangleleft F_{\tau} \tilde{\mathbf{C}}_{nn}^{k} F_{s} \triangleright $$

$$ \mathbf{K}_{L}^{ij\tau s} \quad \text{cost} \ I_{l}^{ij} \triangleleft \left( \sqrt{\frac{b}{b+x}} \frac{x}{F_{\tau,x}(x, z_{\text{top}})} \mathbf{I}_{L} F_{s}(x, z_{\text{top}}) \right) \triangleright x \quad (24) $$

$$ \mathbf{K}_{L\Omega}^{ij\tau s} \quad \text{cost} \ \Omega^2 I_{ly2}^{ij} \triangleleft \left( \sqrt{\frac{b}{b+x}} \frac{x}{F_{\tau,x}(x, z_{\text{top}})} \mathbf{I}_{L} F_{s}(x, z_{\text{top}}) \right) \triangleright x $$

$$ \mathbf{F}_{\Omega}^{\tau s} \quad I_{ly}^{i} \triangleleft F_{\tau} \rho \mathbf{r} \triangleright \Omega^{T} \Omega $$

where

$$ \mathbf{I}_{Ay} \begin{bmatrix} 0 & 0 & 1 \\ 1 & 0 & 0 \\ 0 & 1 & 0 \end{bmatrix} \triangleleft \dots \triangleright \int_{A} \dots dA \triangleleft \dots \triangleright x \int_{b} \dots dx $$

$$ (I_{l}^{ij}, I_{l}^{j,y}, I_{l}^{y,j}, I_{l}^{i,j,y}) \quad \int_{l} (N_{i} N_{j}, N_{i} N_{j,y}, N_{i,y} N_{j}, N_{i,y} N_{j,y}) dy $$

$$ (I_{ly}^{ij}, I_{ly2}^{ij}, I_{\sigma_0}^{ij,y}) \quad \int_{l} (y N_{i} N_{j}, y^{2} N_{i} N_{j}, \sigma_{0} N_{i,y} N_{j,y}) dy $$

To obtain the natural frequencies and the normal modes of the rotor, we solve the homogeneous equation:

$$ \mathbf{M} \ddot{\mathbf{q}} + \mathbf{D}_{\text{TOT}} \dot{\mathbf{q}} + \mathbf{K}_{\text{TOT}} \mathbf{q} = \mathbf{0} \quad (25) $$

Assuming a periodic solution

$$ \mathbf{q} \quad \mathbf{q} e^{i\omega t} \quad (26) $$

substituting Eq. (26) and its derivative in Eq. (25), we obtain

$$ \mathbf{q} e^{i\omega t} [\mathbf{K}_{\text{TOT}} + \mathbf{G}_{\text{TOT}} i\omega \quad \mathbf{M}\omega^2] \quad \mathbf{0} \quad (27) $$

The quadratic eigenvalues problem (QEP) of generic order R in Eq. (27) is written in a classical linear system of order 2×R:

$$ \left\{ \begin{array}{l} \mathbf{M} \ddot{\mathbf{q}} + \mathbf{D}_{\text{TOT}} \dot{\mathbf{q}} + \mathbf{K}_{\text{TOT}} \mathbf{q} = \mathbf{0} \\ \dot{\mathbf{q}} + \dot{\mathbf{q}} = \mathbf{0} \end{array} \right. \quad (28) $$

Now, introducing:
<table><caption>**Table 1**
Dimensions [mm] and material properties [GPa], [kg/m³].</caption><thead><tr><th>Case</th><th>c</th><th>t</th><th>L<sub>w</sub></th><th>r<sub>h</sub></th><th>E<sub>11</sub></th><th>E<sub>22</sub></th><th>E<sub>33</sub></th><th>G<sub>12</sub></th><th>G<sub>31</sub></th><th>G<sub>23</sub></th><th>ν<sub>12</sub></th><th>ν<sub>31</sub></th><th>ν<sub>23</sub></th><th>ρ</th></tr></thead><tbody><tr><td>1</td><td>25.4</td><td>1.6002</td><td>444.5</td><td>63.5</td><td>68.9</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>0.3</td><td>-</td><td>-</td><td>2713</td></tr><tr><td>2</td><td>86.36</td><td>21.59</td><td>863.6</td><td>0.0</td><td>0.0827</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>0.202</td><td>-</td><td>-</td><td>190.99</td></tr><tr><td>3</td><td>76.2</td><td>0.804<sup>a</sup></td><td>305.0</td><td>0.0</td><td>98.0</td><td>7.9</td><td>7.9</td><td>5.6</td><td>5.6</td><td>5.6</td><td>0.28</td><td>0.28</td><td>0.28</td><td>1520</td></tr></tbody></table>

<figure><img src="image_3.png" /><figcaption>**Fig. 1.** First 8 frequencies vs rotational speed for cases 1, 2 and 3.</figcaption></figure>

$$ a \begin{bmatrix} q \\ \dot{q} \end{bmatrix} \dot{a} \begin{bmatrix} \dot{q} \\ \ddot{q} \end{bmatrix} \qquad (29) $$

the equations of motion assume the form:

$$ \frac{R}{T} \frac{1}{i\omega} I = 0 \qquad (30) $$

where:

$$ T^1 R \begin{bmatrix} (K_{TOT})^1 (D_{TOT})^1 (K_{TOT})^1 (M) \\ I & 0 \end{bmatrix} \qquad (31) $$

Since the aerodynamic matrices depend on the reduced frequency 'k', the solutions are found through an iterative procedure [55]. The frequencies are introduced by trial and error for any speed in the considered range until Eq. (30) is obtained. The flutter velocity is determined when the real part of one eigenvalue is null.

## 5. Numerical results

The illustrated approach enables to study the dynamics of wings and blades subjected to aerodynamic loads. Therefore, the following section presents results related to vibrational and flutter

<footnote>a total thickness of the six-layer beam.</footnote>
<table><caption>**Table 2**
Flutter velocities for the three configurations of Table 1.</caption><thead><tr><th>Case</th><th>Theory</th><th>Lamination</th><th>Present</th><th>DLM + 1D-CUF</th><th>Anl. [58]</th><th>EXP. [58]</th></tr></thead><tbody><tr><td>1</td><td>TE3</td><td>Isotropic</td><td>134.9</td><td>131.5<sup>(2.5%)</sup></td><td>-</td><td>-</td></tr><tr><td>2</td><td>TE5</td><td>Isotropic</td><td>64.4</td><td>60.2<sup>(6.98%)</sup></td><td>-</td><td>-</td></tr><tr><td rowspan="4">3</td><td>TE4</td><td>[0<sub>2</sub>/90]<sub>s</sub></td><td>23.5</td><td>23.2<sup>(1.29%)</sup></td><td>21.0<sup>(11.9%)</sup></td><td>25<sup>(6.0%)</sup></td></tr><tr><td>TE4</td><td>[45/ - 45/0]<sub>s</sub></td><td>44.2</td><td>40.3<sup>(9.68%)</sup></td><td>39.0<sup>(13.3%)</sup></td><td>&gt;32</td></tr><tr><td>TE4</td><td>[45<sub>2</sub>/0]<sub>s</sub></td><td>27.8</td><td>26.2<sup>(6.11%)</sup></td><td>27.8<sup>(0.0%)</sup></td><td>28<sup>(0.71%)</sup></td></tr><tr><td>TE4</td><td>[30<sub>2</sub>/0]<sub>s</sub></td><td>26.7</td><td>25.9<sup>(3.1%)</sup></td><td>27.8<sup>(3.96%)</sup></td><td>27<sup>(1.11%)</sup></td></tr></tbody></table>

: results not provided.

(%): percentage difference with respect to the present results (Ref. - Pres.) / Ref.

<figure><img src="image_4.png" /><figcaption>**Fig. 2.** Flutter velocity ($\Omega_f$) and frequency ($f_f$) as functions of the inflow ratio. '-[ ]-': Case 1; '-[ ]-': Case 3 [0/0/90]s. TE3.</figcaption></figure>

<table><caption>**Table 3**
Flutter velocities [rpm] and frequencies [Hz] for the case 3 for various laminations. $\lambda_0=0.05$.</caption><thead><tr><th rowspan="2">Case</th><th rowspan="2">Lamination</th><th colspan="2">TE3</th><th colspan="2">TE4</th></tr><tr><th>Ω<sub>f</sub></th><th>f<sub>f</sub></th><th>Ω<sub>F</sub></th><th>f<sub>F</sub></th></tr></thead><tbody><tr><td rowspan="4">3</td><td>[0<sub>2</sub>/90]<sub>s</sub></td><td>867.45</td><td>26.72</td><td>867.26</td><td>26.72</td></tr><tr><td>[45/ - 45/0]<sub>s</sub></td><td>1621.47</td><td>45.70</td><td>1621.47</td><td>45.67</td></tr><tr><td>[45<sub>2</sub>/0]<sub>s</sub></td><td>1149.73</td><td>32.75</td><td>1124.90</td><td>32.02</td></tr><tr><td>[30<sub>2</sub>/0]<sub>s</sub></td><td>1159.28</td><td>33.23</td><td>1150.69</td><td>32.99</td></tr></tbody></table>

analysis performed on fixed and rotary configurations made of iso-tropic and composite materials.

### 5.1. Rectangular cross section beams

We firstly consider three beams with rectangular cross sections with the geometrical and mechanical properties listed in Table 1.

Cases 1 and 2 are isotropic beams with different material properties (aluminum and foam), whereas the case 3 is a laminated plate with six layers of the equal thickness. Fig. 1 shows the variations of the natural frequencies with the rotational speed for all three cases.

In Fig. 1(a) and (b), the results are compared with the 3D finite element solutions reported in [56]. For the metallic beam (case 1), the TE3 results perfectly match the reference solution, whereas, for the foam beam (case 2), the order of the theory has been increased for describing the significant 3D effects due to the soft material. Fig. 1(c) and (d) show the frequency variations for the case 3 for two symmetric stacking sequences [45<sub>2</sub>/0]<sub>s</sub> and [45/ 45/0]<sub>s</sub>, respectively. It should be noted that, for the first lamination

scheme, the TE3 theory predicts higher values of the 5th and 7th frequencies with respect to the TE4 solution. The corresponding modal shapes are dominated by torsional deformations. Conversely, regarding the [45/ 45/0]<sub>s</sub> scheme, the two Taylor like expansions provide almost equivalent results for all modal shapes.

<figure><img src="image_5.png" /><figcaption>**Fig. 3.** Airfoil NACA 0012. The dimensions are in [mm].</figcaption></figure>
<table><caption>**Table 4**
Material properties of airfoil NACA 0012. Elastic moduli in [GPa] and densities in [kg/m³].</caption><thead><tr><th>Material</th><th>E<sub>11</sub></th><th>E<sub>11</sub></th><th>E<sub>11</sub></th><th>G<sub>12</sub></th><th>G<sub>31</sub></th><th>G<sub>23</sub></th><th>ν<sub>12</sub></th><th>ν<sub>31</sub></th><th>ν<sub>23</sub></th><th>ρ</th></tr></thead><tbody><tr><td>1</td><td>139.0</td><td>11.0</td><td>11.0</td><td>6.05</td><td>6.05</td><td>3.78</td><td>0.313</td><td>0.313</td><td>0.313</td><td>1578</td></tr><tr><td>2</td><td>35.756</td><td>15.427</td><td>15.427</td><td>3.902</td><td>3.902</td><td>3.902</td><td>0.27</td><td>0.27</td><td>0.27</td><td>1826.9</td></tr><tr><td>3</td><td>1.446</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>0.49</td><td>-</td><td>-</td><td>13683.</td></tr><tr><td>4</td><td>0.0827</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>0.2</td><td>-</td><td>-</td><td>190.99</td></tr><tr><td>5</td><td>0.03105</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>0.45</td><td>-</td><td>-</td><td>95.48</td></tr></tbody></table>

-: the material is isotropic.

<figure><img src="image_6.png" /><figcaption>**Fig. 4.** Frequency vs rotational speed for the multi-section NACA 0012 blade. $\theta$ 0.</figcaption></figure>

The flutter analyses have been performed to determine the instability conditions of the three fixed wing configurations. Table 2 lists the flutter velocities computed with the unsteady Theodorsen theory within the CUF framework. For cases 1 and 2, the results are compared with those obtained using the DLM combined with the 1D CUF elements [57]. For case 3, further reference solutions have been taken from [58], where Hollowell et al. studied the aeroelastic behavior of orthotropic plates through analytical and experimental simulations.

The comparisons reveal that the unsteady strip theory predicts higher flutter velocities than those of DLM for the first two cases with a maximum percentage difference equal to about 7%. For case 3, the maximum discrepancy increases up to about 10% for the lamination scheme [45/ 45/0]s. Except for this stacking sequence, it should be noted that the CUF results are close to the experimental data.

In addition, further flutter analyses have been performed with the contributions of the rotational speed. As pointed out in Section 3.1, the distribution of the aerodynamic loads depends on the value of the spacing factor, $\hat{h}$. Therefore, the first analysis has been devoted to evaluating the variations of the instability conditions by varying the inflow ratio, $\lambda_0$, for the metallic and laminated blade (with $Q=4$). The results are shown in Fig. 2.

The x axis reports the logarithmic value of the inflow ratios starting from 0.0242, which corresponds to $\hat{h}=1.0$ for the metallic structure. For the case 1, it is interesting to note that the flutter speed suddenly increases, whereas the frequencies decreases, when $\lambda_0$ overcomes the values 0.03 ($\hat{h} \approx 1.238$) and 0.07 ($\hat{h} \approx 2.89$). Before and after these thresholds, the curves are nearly flat. Conversely, the flutter conditions for the case 2 show a quasi linear dependence with the inflow ratio before the reaching of the asymptotic values. For both structures, for high values of $\lambda_0$, the Loewy approach coincides with the Theodorsen theory. Table 3 reports the instability conditions for the different lamination sequences of the laminated plate assuming the inflow ratio $\lambda_0 = 0.05$.

## 5.2. Rotating multi section blade: NACA 0012

The airfoil of the blade is the NACA 0012, and the chord and length are 0.3454 and 5.186 m, respectively. Fig. 3 shows the blade profile that is constituted by different sections whose material properties are listed in Table 4. The structure is cantilevered and

<table><caption>**Table 5**
Flutter velocities [rpm] and frequencies [Hz] for the NACA 0012.</caption><thead><tr><th rowspan="2"></th><th rowspan="2"></th><th rowspan="2"></th><th colspan="7">θ</th></tr><tr><th>[-45]</th><th>[-30]</th><th>[-15]</th><th>[0]</th><th>[15]</th><th>[30]</th><th>[45]</th></tr></thead><tbody><tr><td rowspan="4">C(k)=1</td><td rowspan="2">TE3</td><td>Ω<sub>f</sub></td><td>424.37</td><td>501.81</td><td>428.66</td><td>212.37</td><td>232.90</td><td>307.86</td><td>305.95</td></tr><tr><td>f<sub>f</sub></td><td>16.03</td><td>16.15</td><td>13.10</td><td>13.50</td><td>17.76</td><td>21.05</td><td>19.14</td></tr><tr><td rowspan="2">TE4</td><td>Ω<sub>f</sub></td><td>378.53</td><td>470.20</td><td>409.85</td><td>212.47</td><td>224.21</td><td>292.59</td><td>289.53</td></tr><tr><td>f<sub>f</sub></td><td>13.69</td><td>15.28</td><td>12.64</td><td>13.00</td><td>17.09</td><td>19.96</td><td>18.30</td></tr><tr><td rowspan="4">Theodorsen</td><td rowspan="2">TE3</td><td>Ω<sub>f</sub></td><td>446.14</td><td>518.62</td><td>455.69</td><td>289.15</td><td>262.31</td><td>347.59</td><td>351.41</td></tr><tr><td>f<sub>f</sub></td><td>15.26</td><td>16.07</td><td>13.48</td><td>12.07</td><td>16.17</td><td>19.38</td><td>18.00</td></tr><tr><td rowspan="2">TE4</td><td>Ω<sub>f</sub></td><td>385.21</td><td>490.92</td><td>437.16</td><td>278.64</td><td>252.38</td><td>330.59</td><td>332.50</td></tr><tr><td>f<sub>f</sub></td><td>13.05</td><td>15.27</td><td>13.00</td><td>11.59</td><td>15.59</td><td>18.39</td><td>17.21</td></tr><tr><td rowspan="4">Loewy</td><td rowspan="2">TE3</td><td>Ω<sub>f</sub></td><td>431.05</td><td>252.10</td><td>105.04</td><td>281.51</td><td>247.13</td><td>331.36</td><td>346.73</td></tr><tr><td>f<sub>f</sub></td><td>14.95</td><td>0.283</td><td>0.151</td><td>11.57</td><td>14.36</td><td>17.64</td><td>17.16</td></tr><tr><td rowspan="2">TE4</td><td>Ω<sub>f</sub></td><td>379.01</td><td>239.01</td><td>103.13</td><td>253.72</td><td>239.01</td><td>315.89</td><td>327.82</td></tr><tr><td>f<sub>f</sub></td><td>13.00</td><td>0.270</td><td>0.150</td><td>9.871</td><td>13.94</td><td>16.80</td><td>16.38</td></tr></tbody></table>
<figure><img src="image_7.png" /><figcaption>Fig. 5. Damping coefficients vs $\dot{m}$. TE3, $\theta$ -15.</figcaption></figure>

it is modeled with 7 4 node beams. Several problem parameters are evaluated including the lamination angle of the skin ply, the aerodynamic theories and the order of TE models.

Firstly, for the ply angle $\theta$ equal to zero, Fig. 4 shows the natural frequencies as a function of the rotational speed by adopting different displacement theories.

The comparison reveals that the three expansions predict similar results for the mode shapes dominated by the bending whereas, slight discrepancies are evident for the torsional deformations (4th and 6th).

Secondly, the flutter analyses have been carried out by using the aerodynamic model proposed in Section 3.1 with different expressions of the lift deficiency function. Assuming that the four blade rotor is operating at a thrust coefficient $C_T$ 0.014, the corresponding values of inflow ratio '$\lambda_0$' and spacing factors '$h$' are 0.084 and 3. Table 5 shows the flutter conditions as functions of the lamination angle. The used displacement theories are TE3 and TE4 expansions.

It is observed that Loewy and Theodorsen theories lead to different results when the angle is 15 and 30. These discrepancies are due to the possibility to have a negative pitch damping, adopting the Loewy approach, which may imply the single degree of freedom flutter in blade torsion. In fact, from Fig. 5, where the damping coefficients are plotted as functions of $\dot{m}$ for $\theta$ 15, it is noted that the pitch damping reaches high negative values when $\dot{m}$ is approaching 0. Conversely, the coefficient related to the flap motion is positive for all values of $\dot{m}$.

Moreover, it is observed significant differences between the flutter velocities obtained through the TE3 and TE4 theories due to the bending torsion coupling that affects the mode shapes. The percentage differences between the two sets of results vary from 2% to 16% depending on the lamination sequence and the considered aerodynamic theory.

## 6. Conclusion

This work aims to present a simple tool for the aeroelastic analysis of rotating and fixed wings. Unsteady aerodynamic theories are coupled with refined beam theories obtained through Carrera Unified Formulation. Numerical simulations have been performed considering rectangular and multi section profiles constituted by isotropic and composite materials. The effects of the aerodynamic

theories, the lamination sequences and the order of structural models on the flutter conditions have been evaluated. The results have been compared, when possible, with theoretical, experimental and numerical solutions. To conclude, the proposed formulation enables

* to study the dynamics of several rotating blade configurations with a significant accuracy and a lower number of degrees of freedom than the 3D solution;

* to predict, if combined with the unsteady aerodynamic theory, the flutter conditions of isotropic and laminated fixed wings;

* to predict the flutter instability of blades in hover with simplified and complex shaped profiles.

Within this work, Carrera Unified Formulation has shown excellent performances in the study of the dynamics and the aeroelasticity of rotors, therefore paving the way for future research that will include:

* extension to nonlinear structural analysis (geometrical as well as material nonlinearities);

* analysis of complex shaped rotating structures (swept tip blades, twisted blades etc.);

* extension to more sophisticated aerodynamic theories (panel methods, CFD)

## References

[1] Cicala P. Le azioni aerodinamiche sui profili di ala oscillanti in presenza di corrente uniforme. Mem R Accad Torino 1934-1935.

[2] Ellenberger G. Bestimmung der luftkräfte auf einen ebenen tragflügel mit querruder. Z angew Math Mech 1936.

[3] Kussner HG. Zusammenfassender bericht über den instationären auftrieb von flugeln. Luftfahrtforsch 1936.

[4] von Karman T, Burgers JM. General aerodynamic theory of perfect fluids. Aerodynamic theory, vol. II; 1943. Originally published in 1934.

[5] Theodorsen T. General theory of aerodynamic instability and the mechanism of flutter. N.A.C.A. report 496, 1935.

[6] Qin Z, Librescu L. Flutter and divergence aeroelastic characteristics for composite forward swept cantilevered wing. J Aircr 1985;22:1001-7.

[7] Guo SJ, Banerjee JR, Cheung CW. The effect of laminate lay-up on the flutter speed of composite wings. Proc IMechE Part G J Aerosp Eng 2003;217:115-22.

[8] Guo SJ, Cheng W, Degang C. Aeroelastic tailoring of composite wing structures by laminate layup optimization. Tech Notes AIAA J 2006;44:3146-9.
[9] Banerjee JR. A simplified method for the free vibration and flutter analysis of bridge decks. J Sound Vib 2003;260:829-45.

[10] Librescu L, Khdeir A. Aeroelastic divergence of swept-forward composite wings including warping restraint effect. AIAA J 1988;26:1373-7.

[11] Gem F, Librescu L. Aeroelastic tailoring of composite wings exhibiting nonclassical effects and carrying external stores. J Aircr 2000;37:1097-104.

[12] Qin Z, Librescu L. Aeroelastic instability of aircraft wings modelled as anisotropic composite thin-walled beams in incompressible flow. J Fluids Struct 2003;18:43-61.

[13] Cesnik CES, Hodges DH, Patil MJ. Aeroelastic analysis of composite wings. In: Proceedings of the 37th structures, structural dynamics and materials conference, Salt Lake City, Utah; 1996.

[14] Patil MJ. Aeroelastic tailoring of composite box beams. In: Proceedings of the 35th aerospace sciences meeting and exhibit, Reno, Nevada; 1997.

[15] Palacios R, Epureanu BI. An intrinsic description of the nonlinear aeroelasticity of very flexible wings. In: Proceedings of the 52nd AIAA/ASME/ASCE/AHS/ASC structures, structural dynamics and materials conference, Denver, Colorado; 2011.

[16] Bennet RM, Edwards JW. An overview of recent developments in computational aeroelasticity. AIAA Paper-98-2421; 1998.

[17] Huttsell L, Schuster D, Volk J, Giesing J, Love M. Evaluation of computational aeroelasticity codes for loads and flutter. AIAA Paper 2001-0569; 2001.

[18] Livne E. Future of airplane aeroelasticity. J Aircr 2003;40:1066-92.

[19] Demasi L, Livne E. Dynamic aeroelasticity of structurally nonlinear configurations using linear modally reduced aerodynamic generalized forces. AIAA J 2009;47:71-90. http://dx.doi.org/10.2514/1.34797.

[20] Loewy RG. Review of rotary-wing v/stol dynamic and aeroelastic problems. J Am Helicopter Soc 1969;14:3-23.

[21] Ham HD. Helicopter blade flutter. Agard report No. 607 (revision of part III of the Agard Manual on aeroelasticity), vol. 10; 1973.

[22] R. Dat. Aeroelasticity of rotary wing aircraft. Agard Lecture Series No. 63 on Helicopter Aerodynamics and Dynamics, vol. 4; 1973.

[23] Loewy RG. A two-dimensional approach to the unsteady aerodynamics of rotary wings. J Aeronaut Sci 1957;24:81-92.

[24] Greenberg JM. Airfoil in sinusoidal motion in a pulsating stream. NACA TN 1326; 1947.

[25] Friedmann P, Yuan C. Effect of modified aerodynamic strip theories on rotor blade aeroelastic stability. AIAA J 1977;15:932-40.

[26] Yuan C, Friedmann P. Aeroelasticity and structural optimization of composite helicopter rotor blades with swept tips. NASA report; 1995.

[27] Piccione E, Bernardini G, Gennaretti M. A finite element model for structural and aeroelastic analysis of rotary wings with advanced blade geometry. In: Proceedings of the COMSOL conference 2010 Paris; 2010.

[28] Bir GS, Chopra I. Aero mechanical stability of rotorcraft with advanced geometry blades. Math Comput Modell 1994;19:159-91.

[29] Friedmann P, Hodges DH. Recent developments in rotary-wing aeroelasticity. J Aircr 1977;14:1027-41. http://dx.doi.org/10.2514/3.58887.

[30] Friedmann P, Hodges DH. Rotary wing aeroelasticity – a historical perspective. J Aircr 2003;40:1019-46. http://dx.doi.org/10.2514/2.7216.

[31] Friedmann P. Rotary-wing aeroelasticity: current status and future trends. AIAA J 2004;42:1953-72. http://dx.doi.org/10.2514/1.9022.

[32] Conlisk AT. Modern helicopter aerodynamics. Annu Rev Fluid Mech 1997;29:515-67.

[33] Ormiston RA. Rotor dynamic considerations for large wind power generator systems. In: Proc. wind energy conversion systems workshop, national science foundation, NSF/RA/W-73-006; 1973.

[34] Friedmann PP. Aeroelastic modeling of large wind turbines. J Am Helicopter Soc 1976;21:17-27.

[35] Kaza KRV, Hammond CE. Investigation of flap-lag stability of wind turbine rotors in the presence of velocity gradients and helicopters in forward flight. In: 17th Structures, structural dynamics, and materials conference; 1976, May.

[36] Friedmann PP. Aeroelastic stability and response analysis of large horizontal-axis wind turbines. J Ind Aerodyn 1980;5:373-401.

[37] Rasmussen F, Hansen K, Thomsen MH, Larsen TJ, Bertagnolio F, Johansen J, et al. Present status of aeroelasticity of wind turbines. Wind Energy 2003;6:213-28. http://dx.doi.org/10.1002/we.98.

[38] Hansen MOL, Sorensen JN, Voutsinas S, Sorensen N, Madsen HAa. State of the art in wind turbine aerodynamics and aeroelasticity. Progr Aerosp Sci 2006;42:285-330. http://dx.doi.org/10.1016/j.paerosci.2006.10.002.

[39] Carrera E, Filippi M, Zappino E. Laminated beam analysis by polynomial, trigonometric, exponential and zig-zag theories. Eur J Mech A/Solids 2013;41:58-69. http://dx.doi.org/10.1016/j.euromechsol.2013.02.006.

[40] Carrera E, Filippi M, Zappino E. Free vibration analysis of laminated beam by polynomial, trigonometric, exponential and zig-zag theories. J Compos Mater 2013.

[41] Daoud SM, Carrera E, Zenkour AM, Al Khateeb SA, Filippi M. Free vibration of FGM layered beams by various theories and finite elements. Compos Part B 2014;59.

[42] Carrera E, Petrolo M, Zappino E. Performance of cuf approach to analyze the structural behavior of slender bodies. J Struct Eng 2012;138:285-98.

[43] Carrera E, Zappino E, Filippi M. Free vibration analysis of thin-walled cylinders reinforced with longitudinal and transversal stiffeners. J Vib Acoust 2013;135.

[44] Carrera E, Filippi M, Zappino E. Free vibration analysis of rotating composite blades via carrera unified formulation. Compos Struct 2013;106:317-25. http://dx.doi.org/10.1016/j.compstruct.2013.05.055.

[45] Carrera E, Filippi M, Zappino E. Analysis of rotor dynamic by one-dimensional variable kinematic theories. J Eng Gas Turbines Power 2013;135:092501. http://dx.doi.org/10.1115/1.4024381.

[46] Carrera E, Filippi M. Variable kinematic one-dimensional finite elements for the analysis of rotors made of composite materials. ASME J Gas Turb Pwr 2014;136:092501. http://dx.doi.org/10.1115/1.4027192.

[47] Carrera E, Filippi M. Vibration analysis of thin/thick, composites/metallic spinning cylindrical shells by refined beam models. J Vib Acoust 2015;137:031020. http://dx.doi.org/10.1115/1.4029688.

[48] Varello A, Lamberti A, Carrera E. Static aeroelastic response of wing-structures accounting for in-plane cross-section deformation. Int J Aeronaut Space Sci 2013;14.

[49] Carrera E, Varello A, Demasi L. A refined structural model for static aeroelastic response and divergence of metallic and composite wings. CEAS Aeronaut J 2013;4.

[50] Filippi M, Carrera E. Aeroelastic and mechanical hierarchical aeroelastic analysis of composite wings. Mech Adv Mater Struct, in print.

[51] Carrera E, Giunta G, Petrolo M. Beam structures. Classical and advanced theories. Wiley; 2011.

[52] Bisplinghoff RL, Ashley H, Halfman RL. Aeroelasticity. Dover Publications, INC.; 1996.

[53] Rodden WP. Theoretical and computational aeroelasticity. Crest Publishing; 2012.

[54] Bielawa RL. Rotary wing structural dynamics and aeroelasticity. AIAA; 1992.

[55] Nguyen N. Flight dynamics of flexible aircraft with aeroelastic and inertial force interactions. In: AIAA atmospheric flight mechanics conference, vol. 6045; 2009.

[56] Truong KV, Yeo H, Ormiston RA. Structural dynamics modeling of rectangular rotor blades. Aerosp Sci Technol 2013;30:293-305.

[57] Petrolo M. Flutter analysis of composite lifting surfaces by the 1d carrera unified formulation and the doublet lattice method. Compos Struct 2013;95:539-46. http://dx.doi.org/10.1016/j.compstruct.2012.06.021.

[58] Hollowell SJ, Dugundji J. Aeroelastic flutter and divergence of stiffness coupled, graphite/epoxy cantilevered plates. J Aircr 1984;21:69-76. http://dx.doi.org/10.2514/3.48224.