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

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.


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 com bined with simplified structural theories to develop reliable aeroe lastic 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 intrin sic 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 aerody namic theories were proposed in [21,22], respectively. In that per iod, 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 varia tion 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 posi tion; 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 exten sively 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].

The structural model: Carrera Unified Formulation
CUF states that the displacement field, u T ðx; y; z; tÞ ðu x ; u y ; u z Þ T , is an expansion of generic functions, F s ðx; zÞ for the displacement vector, u s ðyÞ: uðx; y; z; tÞ F s ðx; zÞu s ðy; tÞ s 1; 2; . . . ; T where T is the number of terms in the expansion and, according to Einstein's generalized notation, s indicates summation. In this work, Eq. (1) consists of Taylor like expansions denoted as 'TE'. For example, the third order displacement field (TE3) is: 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: p f zz xx xz g T r p fr zz r xx r xz g T n f zy xy yy g T r n fr zy r xy r yy g T ð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: r pCpp p þC pn n r nCnp p þC nn n ð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 arẽ 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 u s ðy; tÞ N i ðyÞq si ðtÞ ð 6Þ where N i ðyÞ are the shape functions and q si ðtÞ is the nodal displace ment vector

2-D steady and unsteady aerodynamic theories for fixed-wing
Despite steady and unsteady strip theories were already com bined 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 a(t). As derived in [52], the final expres sion of the lift for unit of span is where 'U' is the free stream velocity, 'b' the semichord, 'q a ' the air density and 'CðkÞ' the deficiency lift function, which depends on the so called reduced frequency k xU b . The term 'a' defines the posi tion 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Þ FðkÞ þ iGðkÞ H where 'L w ' is the length of the wing, 'c' is the mean chord and, 'c la0 ' is the lift curve slope equal to 2p. Finally, since the quantity 'bp' may be approximated with the integral bþx dx, which reproduces the distribution of pressure over a slightly inclined, thin and uncam bered airfoil, Eqs. (8) and (10) L u 2p AR cosðKÞ

Unsteady aerodynamic theory for blade
According to [54], the general expression for the lift distribution on a rotary wing yields dL dy where the reduced frequency 'k' is 2xc 3XR and, 'R' is rotor radius. Although this formulation is limited as it relies on Theodorsen prob lem, 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 increas ing 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, where 'k 0 ' C T =2 p ' (where C T is the thrust coefficient) is the inflow ratio and, 'r' is the generic radial distance.

Equations of motion
The general form of Hamilton's Principle states that the varia tion 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 where d refers to kinematically admissible perturbations during the indicated time interval and T: total kinetic energy of the system U + U r0 : total potential energy of the system, including both strain energy ('U') and potential of any conservative external forces ('U r0 ') L 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 where 'q' 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 The external work generated by the lift can be written as where 'z top ' 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 ext cost The subscripts 'L' and 'X' denote the contributions due to the aero dynamic load and rotation, respectively. The compact expressions of the fundamental nuclei yield ; I i; y j l ; I i; y j; y l Z l N i N j ; N i N j; y ; N i; y N j ; N i; y N j; y dy I i j ly ; I i j l y 2 I i; y j; y lr 0 Z l y N i N j ; y 2 N i N j ; r 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: Assuming a periodic solution substituting Eq. (26) and its derivative in Eq. (25), we obtain The quadratic eigenvalues problem (QEP) of generic order R in Eq. (27) is written in a classical linear system of order 2ÂR: Now, introducing: the equations of motion assume the form: where: 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.

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  analysis performed on fixed and rotary configurations made of iso tropic and composite materials.

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 prop erties (aluminum and foam), whereas the case 3 is a laminated plate with six layers of the equal thickness. Fig. 1 shows the varia tions 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)     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 com bined 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 experimen tal 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, 'ĥ'. Therefore, the first analysis has been devoted to evaluating the variations of the instability condi tions by varying the inflow ratio, 'k 0 ', for the metallic and lami nated 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ĥ 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 k 0 overcomes the values 0.03 (ĥ ' 1:238) and 0.07 (ĥ ' 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 k 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 k 0 0:05.

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   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 h equal to zero, Fig. 4 shows the natural frequencies as a function of the rotational speed by adopting differ ent displacement theories.
The comparison reveals that the three expansions predict simi lar 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 expres sions of the lift deficiency function. Assuming that the four blade rotor is operating at a thrust coefficient C T 0:014, the corre sponding values of inflow ratio 'k 0 ' and spacing factors 'ĥ' 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 dif ferent results when the angle is 15 and 30. These discrepancies are due to the possibility to have a negative pitch damping, adopt ing the Loewy approach, which may imply the single degree of freedom flutter in blade torsion. In fact, from 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.

Conclusion
This work aims to present a simple tool for the aeroelastic anal ysis 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, experimen tal and numerical solutions. To conclude, the proposed formulation enables to study the dynamics of several rotating blade configura tions with a significant accuracy and a lower number of degrees of freedom than the 3D solution; to predict, if combined with the unsteady aerodynamic the ory, the flutter conditions of isotropic and laminated fixed wings; to predict the flutter instability of blades in hover with sim plified and complex shaped profiles.
Within this work, Carrera Unified Formulation has shown excel lent performances in the study of the dynamics and the aeroelas ticity 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)