New Method for Analyzing the Flutter Stability of Hingeless Blades with Advanced Geometric Configurations in Hovering

A new method used to analyze the aeroelastic stability of a helicopter hingeless blade in hovering has been developed, which is especially suitable for a blade with advanced geometric configuration. This method uses a modified doublet-lattice method (MDLM) and a 3-D finite element (FE) model for building the aeroelastic equation of a blade in hovering. Thereafter, the flutter solution of the equation is calculated by the V-g method, assuming blade motions to be small perturbations about the steady equilibrium deflection. The MDLM, which is suitable to calculate the unsteady aerodynamic force of nonplanar rotor blade in hovering, is developed from the doublet-lattice method (DLM). The structural analysis tool is the commercial software ANSYS. The comparisons of the obtained results against those in the literatures show the capabilities of the MDLM and the method of structural analysis. The flutter stabilities of swept tip blades with different aspect ratios are analyzed using the new method developed in this work and the usual method on the basis of the unsteady strip theory and beam model. It shows that considerable differences appear in the flutter rotational velocities with the decrease of the aspect ratio. The flutter rotational velocities obtained by the present method are evidently lower than those obtained by the usual method.


Introduction
Hovering is the most important flight state of a helicopter. The flutter stability of hovering is directly related to flight safety. Numerous methods on stability characteristics of hingeless rotors in hovering have been developed to obtain an accurate flutter stability boundary. The aeroelastic equations of their methods are usually obtained by aerodynamic models, which are formed by quasisteady [1] or unsteady strip theories [2], and the structural dynamic models, which are formed by moderate [3] and large deflection beam [4] theories or other refined beam theories [5][6][7], considering the geometrical nonlinearity. Generally, for the flutter solution of the aeroelastic equations, the nonlinear steady state is obtained first and then the flutter solution, assum-ing that the flutter motion is a small perturbation about the equilibrium position [1][2][3][4][5][6][7][8]. The blades of aforementioned researches are made of isotropic [3] or composite [8] material with high aspect ratios and simple planform shapes. With the development of computational fluid dynamics (CFD), the method of rotorcraft aeroelastic tight coupling of computational structural dynamics (CSD) and CFD has been used for rotor blade aeroelastic stability analysis recently [9,10]. Although the CFD/CSD coupled method can get better predictions of rotor stability, the computational cost of this method is very high [11], and CSD models are mostly built with beams mentioned above.
The majority of manufacturers choose blades with advanced geometric configurations, such as the ERATO [12] blade with reduced Blade Vortex Interaction (BVI) noise signature to increase the flight speed and have good aeroacoustic behavior. The advanced blade with a complicated tip and a small aspect ratio deviates largely from the classical rectangular one. The modal analysis using the beam model of blades may be unsuitable due to the deviation of a planform shape. Truong [13] studied the dynamic characteristics of the ERATO using 3-D finite element and 1-D beam analyses in MSC software. Natural frequencies are calculated at various rotor rotational speeds. The predictions of the 3-D finite element model that correlate with the experimental results are more precise than those of the (1-D + 2-D) model, within a 5% relative error. Yeo [14] et al. are motivated by the above research to study the difference between 1-D beam and 3-D finite element analyses for advanced geometry blades, which have tip sweep, tip taper, and planform variations near the root and physics behind them. Differences occur when coupling generated from geometry (tip sweep) or material (composite) emerges, especially for highfrequency modes. Thereafter, Kee and Shin [15] use an eighteen-node solid-shell finite element to model the blade structures and give the differences of natural frequencies between the one-dimensional beam and threedimensional solid models. As the blade aspect ratio decreases, considerable differences appear in the bending and torsion modes. Evidently, the accuracy of modal calculation directly determines the reliability of flutter analysis. Here comes the conclusion that the commonly used aeroelastic model constructed by strip theory and the 1-D beam model is unsuitable for the flutter stability analysis of blades with advanced geometric configurations in hovering. It is also easy to conclude that the CFD/CSD tight coupling approach using 3-D solid structure models instead of beam models will increase computing costs inevitably when there is already a large amount of computation. A new high-efficiency approach of advanced blade flutter analysis is sorely necessary. So far, no relevant literature has been reported on that.
A new method of rotor flutter analysis, using a 3-D finite element model and a modified doublet-lattice method (MDLM) to establish the aeroelastic equation, is developed to solve the abovementioned problems. ANSYS is applied to the numerical analysis of blade structure, and the blades are discretized by the 3-D solid element. In aerodynamic calculation, compared with the aerodynamic models mentioned above, one of the major advantages of the DLM [16] is that motions of parts of the wing section can be considered, which makes the DLM more suitable for being coupled with the 3-D finite element model to establish the aeroelastic model than the strip theory. The other one is that the DLM is much more effective than the CFD method. However, the DLM cannot be used for rotor flutter analysis directly because the inflows change over the location of the blade and the pretwist cannot be ignored. The MDLM, which overcomes the above difficulties to calculate the unsteady aerodynamic of rotor blades in hovering, is developed from the DLM. The method of solving the aeroelastic equations is basically the same as the abovementioned one. The nonlinear steady state is obtained first to linearize the flutter equations, based on the small perturbation assumption. Thereafter, the linearized flutter equations are solved through the V-g method [2]. The structural analysis method and MDLM of rotor blades are verified by literature examples. Thereafter, the new and usual methods of rotor flutter analysis are used for solving the flutter rotational velocities of advanced blades in hovering, which have different ratio aspects. The comparison of these two methods is conducted systematically to understand the differences. The DLM assumes the lifting surfaces as a planar plate without pretwist angles; that is, in the aerodynamic coordinate frame of the sending panel, the chordwise component of the normal vector of the receiving panel is equal to zero [17]. Pretwist angles are widely used in helicopter rotor blades, which have a great influence on aerodynamic forces. Hence, the DLM's assumption is unsuitable for rotor blades. On the basis of the DLM, the MDLM is developed for calculating unsteady aerodynamic loads on nonplanar lifting surfaces in subsonic flows.

Analysis
Normal wash velocity [18] of receiving point ðx, y, zÞ induced by the surface motion can be written as follows: and M ∞ are the pressure loading, velocity pressure number, freeflow velocity, and Mach number at sending point ð b ξ, b ζ, b ηÞ, respectively. ω is the circular frequency of vibration of the blade. The kernel function [18] can be expressed as follows: where Partial differentiations ∂/∂n c and ∂/∂n d are taking the derivative of the receiving point and the sending point. ∂

2
International Journal of Aerospace Engineering /∂n c and ∂/∂n d can be expressed by the following formulations: In the case of small perturbations, the following relations are obtained: cos ðn c , xÞ cannot be ignored due to blade elastic torsion and pretwist. Formulation (4) can be rewritten as follows: The kernel function K can be also written as follows: where Because and using ψ s for representing the anhedral of the receiving point and ψ d for representing the anhedral of the sending point, we can write certain equations as follows: And then According to the following equations: and substituting Equations (11) and (12) into Equation (7), the kernel function can be rewritten as follows: 3 International Journal of Aerospace Engineering The computational procedure of K 1 and K 2 are given by Landahl [18]. ∂K 1 /∂xof Equation (13) can be calculated by taking the partial derivative of K 1 with respect to x. The results are given as follows: where The constants a n and c 0 of Equation (16) are given by Kalman et al. [19].
If u 1 ≥ 0, then I 1 , I 2 , and ∂I 1 /∂x will be written as follows: − nc 0 a n ∂u 1 ∂x a n e −nc 0 u 1 International Journal of Aerospace Engineering And if u 1 < 0, then I 1 , I 2 , and ∂I 1 /∂x can be rewritten as follows:  (1) and (13), the normal wash velocity w ij of the control point can be given as follows: In Equation (19), U j is the freestream velocity of the jth panel, S j is the area of the jth panel, and Δ p j is the pressure difference per unit area of the jth panel. The dynamic pressure is given by q dj = 0:5ρ ∞ U 2 j , where ρ ∞ is the air density. Let the pressure difference on the unit length of 1/4 chord of the jth panel be Δ p. To make the distribution of Δ p on 1/4 chord equivalent to that of Δ p j on the jth panel, their relationship is as follows: where l j is the length of the 1/4 chord of the jth panel. Here, S j = ðl j cos χ j ÞΔx j , where χ j is the sweep angle of the doublet line and Δx j is the average chord. The pressure difference across surface Δ p j can be written as Substituting S j and Δ p j into Equation (20), Δ p can be rewritten as follows: Transform the area integral of Equation (19) into a line integral, and then the normal wash velocity is obtained as follows: where The downwash factor D ij can be obtained by the steady downwash factor D 0ij and incremental oscillatory downwash factors D 1ij and D 2ij , where The steady downwash factor D 0ij is derived from horseshoe vortex considerations. The difference from the traditional method given by Hedman [20] is that the infinite free vortex is replaced by the wake shed from the blade. The induced velocities of the wake on the panels are calculated by the classical blade element momentum theory and considered to be steady and uniform along the blade length [21]. For brevity, the derivation processes of incremental oscillatory downwash factors D 1ij and D 2ij are not presented, which can be referred to the literature [20]. The numerical form of Equation (1) in matrix notation becomes The differential pressure coefficient vector ½Δc p can be obtained by ½Δc p = ½D −1 ½w. For the jth panel, the pressure difference located at the 1/4 chord point at the midspan of the panel can be obtained byF aj = q dj ΔS j Δc pj . Thereafter, the aerodynamic load vector is obtained in the form of a matrix by the following: where where ½q s is the displacement vector and ½F s is the structural load vector. M ss and K ss are the structural mass 5 International Journal of Aerospace Engineering matrix and stiffness matrix, respectively, and the functions of rotational speed Ω and q s . Equation (28) is transformed into the modal space by writing where Φ s is the modal matrix and ½p s is the vector of eight generalized coordinates in the modal space. Substituting Equation (29) into Equation (28), the normal equation can be rewritten as follows: where 3. Aeroelastic Modeling. Given the difference of meshing between the aerodynamic and the structure models, an interconnecting approach must be established between them for aeroelastic analysis, and it can be made using a spline matrix G as [22]. The transformation from the structural grid node motions ½q s into the aerodynamic grid node motions ½w can be given by the following: The transformation formula of the aerodynamic ½F a and the equivalent value acting on a structural grid ½F s can be obtained by the following: Substituting Equations (31) and (32) into Equation (26), we can obtain the following: where Q is a function of k and rotor rotational velocity Ω and Combining Equations (28), (33), and (34), the aeroelastic equation is obtained by Transforming Eq. (35) to the modal space, the equation can be rewritten as follows: where Qðk, ΩÞ = Φ T s Qðk, ΩÞΦ s .

Nonlinear Steady-State Solution.
The conclusion of the literature [3] shows that the static deformation of the blade will change the structural stiffness, due to the consideration of geometrical nonlinearity, and the aerodynamic force distribution, which will affect the blade flutter stability. Hence, the first step in solving Equation (35) is to get the steady trim solution in hovering. Dropping the acceleration term, we can write Equation (35) as follows: where stiffness matrix K ss 0 is a function of the steady displacement ½q s0 and Ω and structural load vector ½F s0 is given by the following: The velocity vector ½w 0 of Equation (38) is a vector consisting of the normal downwash velocity, which is the normal component of a combined velocity of freestream and induced inflow of the wake, at each panel. For solving the nonlinear aeroelastic equation (37), an iterative algorithm is established in the Isight platform integrating the steady aerodynamic calculation part of self-made MDLM program and the ANSYS software. The process of the algorithm is as follows: (1) The initial displacements of the structural grid nodes are set by ½q s00 = 0. The displacements of the structural grid nodes of the upper wing surface are taken to represent the overall displacement of the blade. The V-g method is used to calculate the blade flutter rotational velocity Ω F and flutter frequency (2) In the i th iteration, where i = 1, 2, 3, ⋯, n, the displacements of the structural grid nodes are ½q s0i−1 .
The location of structural and aerodynamic grid nodes is refreshed. The normal downwash velocity vector ½w 0 and the steady aerodynamic influence coefficient matrix ½D 0 −1 are recalculated due to the change of the normal direction and the coordinates of sending and receiving points at each panel (3) The interpolation matrix G as is obtained using the Thin Plate Spline (TPS) interpolation method [23]. The steady structural load ½F s0 is obtained solving Equation (38) (4) The structural grid node displacement vector ½q s0i is calculated using the static structural analysis of ANSYS software when the structural load is applied to the upper wing surface (5) If j½q s0i − ½q s0i−1 j < ε, then, the computation converges, and if it does not, then, Steps (2)-(4) are repeated After the calculation converges, the structure displacement ½q s0 , the location of structural and aerodynamic grid nodes, and the structural load ½F s0 in a static trim state are obtained. Given the large deflection and the geometric nonlinearity, the natural modes of vibration Φ s0 , the mass matrix M ss 0 ðΩ, q s0 Þ, and the stiffness matrix K ss 0 ðΩ, q s0 Þ about the equilibrium position are obtained by the prestress modal analysis of ANSYS. 6 International Journal of Aerospace Engineering 2.5. Flutter Solution. The flutter motion ½q s is assumed to be a small perturbation ½q s about the equilibrium position ½q s0 and can be written as follows: Substituting Equation (39) into Equation (35), we can linearize the flutter equation, based on the assumption of small perturbation, by the following: Simplified by Equation (37), Equation (40) can be rewritten as follows: Transforming Equation (41) into the modal space, the equation can be rewritten as follows: where , and ½p s = Φ T s0 ½q s . Based on the assumption of harmonic motion (½p s = ½p ss e iωt ) transforming Equation (42) into a frequency domain, the aeroelastic equation can be rewritten as follows: Equation (43) is solved in the frequency domain by the V-g method to obtain the flutter rotational velocity of the rotor blade in hovering. A specific process can refer to Ref. 24.
Given the effect of centripetal force, the mass and stiffness matrixes of the blade are the functions of rotational velocity. To match the current rotational velocity to the flutter rotational velocity, we can use an iterative algorithm as follows: (1) In the i th iteration, where i = 1, 2, 3, ⋯, n, for the rotational velocity Ω i and the collective pitch θ 0 , the static steady-state solution is carried out, and the specific steps can be referred to the abovementioned After the calculation converges, the flutter rotational velocity Ω F and flutter frequency f F is obtained for the current blade pitch.

MDLM Validation.
In this work, the steady-state deformation and flutter stability of a hovering rectangular hingeless blade are analyzed. The aeroelastic model is established by combining the moderate deflection beam theory with the MDLM.
The blade is made of isotropic material. Nondimensional natural frequencies in the flap, lead-lag, and torsion directions are ω w = 1:15, ω v = 1:5, and ω ϕ = 3:0 or 5.0, respectively. Table 1 shows other blade nondimensional parameters [3]. Figure 1 exhibits the steady tip deflection w tip , flap and lead-lag bending deflection v tip , and elastic twist angle θ tip with respect to collective pitch θ 0 and the results of Sivaneri and Chopra [3]. The agreement between the two results is good. Figure 2 shows the flutter stability boundary respect with β p and the results of Sivaneri and Chopra [3]. The flutter solution is calculated by the V-g method. A general agreement emerges between the two results.

Rotational Modal Analysis.
In the present study, the natural frequencies of the swept tip blades were investigated in ANSYS, and the results were compared with those of the experiment [25] to assess the availability of ANSYS for analyzing the dynamic characteristics of rotating hingeless isotropic rotor blades. The blade parameters are the same as those of the beams tested by Epps and Chandra [25]. The radius of the rotor is 1,016 mm, and the hub length is 63.5 mm. In addition, the length of the swept tip segment is 152.5 mm, the width is 25.4 mm, and the thickness is 1.6 mm. The tip sweep angles Λ are 15 and 45 degrees. The isotropic material of the blades is aluminum. The model is 3-D FE model of 2908 nodes and of 360 elements, and the meshing is shown in Figure 3. Figures 4 and 5 show the natural frequencies with respect to the rotating speeds. The results show that the present predictions have a good correlation with the results of the experiment.

Flutter Solution.
The flutter rotational velocities of hingeless rotor blades that have swept tips and different aspect ratios are obtained for different collective pitches by the present method mentioned above and the usual   [4]. The airfoil section of the blade is naca8h12 and is shown in Figure 6. In Figure 6, the elastic axis position is ea, and the distance to the leading edge of the ea is XE. Table 2 shows the parameters of the section. There exists a destabilizing offset, XA, between the aerodynamic center and the elastic axis.
The isotropic material is aluminum. Figure 7 shows the blade platform, which is based on the blade tested by Harder and Desmarais [22] and includes a hub. The radius of the rotor blade is R = 4:252 m. Two hub lengths (hi = 0:367 m and 1.5645 m) are considered in investigating the effect of aspect ratios (AR = 18:07 and 12.5). The length of the swept tip segment is n = 0:6378 m. The tip sweep angle Λ is 45 degrees. Table 3 shows natural frequencies of the blade (AR = 18:07) at the rotational velocity 400 rpm, as the total number of elements is varied from 4880 to 6710. The results show, for the case considered, that 6100 elements are sufficient for reasonably good convergence. Figures 8-9 show the flutter rotational velocities, and Figures 10-11 show the V-g diagrams, with respect to various aspect ratios, where Ω = 500 rpm. The results show that the flutter rotational velocities obtained by the present method are evidently lower than those obtained by the usual method. The differences in flutter rotational  9 International Journal of Aerospace Engineering velocities between the two methods increase with the decrease of the aspect ratio. When AR = 18:07, the results of the two methods are close, and the differences between the two methods are less than 10%. When AR = 12:5, the differences between them are 10% to 20%.
The cause of the differences may mainly consist of the following two parts: (1) The flutter of both conditions occurs between the 2 nd mode and the 3 rd mode. The modal analysis of these two modes directly affects the flutter result. Figures 12-13 show the frequencies of these modes about steady deflected positions with respect to the rotating speeds (theta 0 = 6°). Figures 14-15 show the shapes of these modes of the blade (AR = 12:5). The 2 nd mode is the flap and torsion coupled mode, and the 3 rd mode is the lead-lag and torsion coupled mode. The differences in modal analysis between the 1-D beam model and the 3-D finite element model, increasing with the decrease of the aspect ratio, because of strong coupling, generated from geometry, between flap, lag and torsion [12], may cause differences of flutter results