Abstract
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.
1. 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–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, assuming that the flutter motion is a small perturbation about the equilibrium position [1–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 () 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 high-frequency 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 three-dimensional 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.
2. Analysis
2.1. Modified Doublet Lattice Method (MDLM)
2.1.1. Correction of the Nonplanar Hypothesis
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.
Normal wash velocity [18] of receiving point induced by the surface motion can be written as follows: where , , , , and are the pressure loading, velocity pressure number, free-flow velocity, and Mach number at sending point , respectively. is the circular frequency of vibration of the blade. The kernel function [18] can be expressed as follows: where
Partial differentiations and are taking the derivative of the receiving point and the sending point. and can be expressed by the following formulations:
In the case of small perturbations, the following relations are obtained:
cannot be ignored due to blade elastic torsion and pretwist. Formulation (4) can be rewritten as follows:
The kernel function can be also written as follows: where
Because and using for representing the anhedral of the receiving point and 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: where
The computational procedure of and are given by Landahl [18]. of Equation (13) can be calculated by taking the partial derivative of with respect to . The results are given as follows: where
The constants and of Equation (16) are given by Kalman et al. [19].
If , then , , and will be written as follows:
And if , then , , and can be rewritten as follows:
2.1.2. Application of the MDLM on Blade
In the MDLM, the lifting surface of a blade is divided into trapezoidal panels, which have different freestream velocities. The freestream velocity and the reduction frequency are assumed to be constant within the same panel to satisfy the basic assumption and condition of the MDLM based on the DLM. Each panel contains a line of constant strength point doublets along the quarter-chord line and a control point located at the 3/4 chord point at midspan. Substituting the coordination of the control point and the doublet point into Equations (1) and (13), the normal wash velocity of the control point can be given as follows:
In Equation (19), is the freestream velocity of the th panel, is the area of the th panel, and is the pressure difference per unit area of the th panel. The dynamic pressure is given by , where is the air density. Let the pressure difference on the unit length of 1/4 chord of the th panel be . To make the distribution of on 1/4 chord equivalent to that of on the th panel, their relationship is as follows: where is the length of the 1/4 chord of the th panel. Here, , where is the sweep angle of the doublet line and is the average chord. The pressure difference across surface can be written as , where is the coefficient of the pressure difference. Substituting and into Equation (20), 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 can be obtained by the steady downwash factor and incremental oscillatory downwash factors and , where
The steady downwash factor 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 and 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 can be obtained by . For the th panel, the pressure difference located at the 1/4 chord point at the midspan of the panel can be obtained by. Thereafter, the aerodynamic load vector is obtained in the form of a matrix by the following: where
2.2. Structure Modeling
Considering the geometric nonlinearity, the nonlinear forced undamped vibration equation of the blade can be simply obtained by where is the displacement vector and is the structural load vector. and are the structural mass matrix and stiffness matrix, respectively, and the functions of rotational speed and . Equation (28) is transformed into the modal space by writing where is the modal matrix and 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 and .
2.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 [22]. The transformation from the structural grid node motions into the aerodynamic grid node motions can be given by the following:
The transformation formula of the aerodynamic and the equivalent value acting on a structural grid can be obtained by the following:
Substituting Equations (31) and (32) into Equation (26), we can obtain the following: where is a function of 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 .
2.4. 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 is a function of the steady displacement and and structural load vector is given by the following:
The velocity vector 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 . 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 and flutter frequency(2)In the th iteration, where , the displacements of the structural grid nodes are . The location of structural and aerodynamic grid nodes is refreshed. The normal downwash velocity vector and the steady aerodynamic influence coefficient matrix 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 is obtained using the Thin Plate Spline (TPS) interpolation method [23]. The steady structural load is obtained solving Equation (38)(4)The structural grid node displacement vector is calculated using the static structural analysis of ANSYS software when the structural load is applied to the upper wing surface(5)If , then, the computation converges, and if it does not, then, Steps (2)–(4) are repeated
After the calculation converges, the structure displacement , the location of structural and aerodynamic grid nodes, and the structural load in a static trim state are obtained. Given the large deflection and the geometric nonlinearity, the natural modes of vibration , the mass matrix , and the stiffness matrix about the equilibrium position are obtained by the prestress modal analysis of ANSYS.
2.5. Flutter Solution
The flutter motion is assumed to be a small perturbation about the equilibrium position 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 .
Based on the assumption of harmonic motion () 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 th iteration, where , for the rotational velocity and the collective pitch , the static steady-state solution is carried out, and the specific steps can be referred to the abovementioned(2)The V-g method is used to calculate the blade flutter rotational velocity and flutter frequency (3)If , the computation converges, and if it does not, and Steps (1)–(2) are repeated
After the calculation converges, the flutter rotational velocity and flutter frequency is obtained for the current blade pitch.
3. Presentation of Results
3.1. 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 , , and or 5.0, respectively. Table 1 shows other blade nondimensional parameters [3].
Figure 1 exhibits the steady tip deflection , flap and lead-lag bending deflection , and elastic twist angle with respect to collective pitch and the results of Sivaneri and Chopra [3]. The agreement between the two results is good.
(a)
(b)
(c)
Figure 2 shows the flutter stability boundary respect with 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.
3.2. 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.
3.3. 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 method constructed by unsteady strip aerodynamic theory and large deflection beam theory [4]. The airfoil section of the blade is naca8h12 and is shown in Figure 6. In Figure 6, the elastic axis position is , and the distance to the leading edge of the 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 m. Two hub lengths ( and 1.5645 m) are considered in investigating the effect of aspect ratios ( and 12.5). The length of the swept tip segment is m. The tip sweep angle is 45 degrees.
Table 3 shows natural frequencies of the blade () 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 . 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 velocities between the two methods increase with the decrease of the aspect ratio. When , the results of the two methods are close, and the differences between the two methods are less than 10%. When , 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 2nd mode and the 3rd 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 (). Figures 14-15 show the shapes of these modes of the blade (). The 2nd mode is the flap and torsion coupled mode, and the 3rd 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(2)In the modal analysis, the blade sections of the 3-D finite element model have deformations, which are not included in the strip theory aerodynamic loadings. But the motions of parts of sections can be considered in the MDLM. Figures 16–17 show the displacements of the 3-D finite element nodes at in the Flap/torsion coupled mode about the steady deflected positions. The results show that the deformations of sections increase with the decrease of aspect ratio
(a) Front view
(b) Top view
4. Conclusion
(1)The results of steady-state deflection and the flutter stability boundaries show that the MDLM is reliable in the steady and unsteady aerodynamic calculations(2)The dynamic characteristics of rotating blades obtained by the modal analysis of ANSYS show good agreement with the experiment(3)The flutter solutions of the blades with advanced geometric configurations are obtained by the new method, based on the MDLM and ANSYS-based structural analysis, and the usual method, based on strip theory and the beam theory model. The flutter rotational velocities obtained by the present method are lower than those obtained by the usual method. The differences of flutter rotational velocities between the two methods increase with the decrease of the aspect ratio. The distinctions may arise mainly from the differences of modal analysis between the 1-D beam model and the 3-D finite element model and the warping of the blade sections. Hence, the method of this paper is more suitable for the flutter stability analysis of advanced geometric blades than the usual method
Nomenclature
| : | Pressure loading |
| : | Velocity pressure number |
| : | Circular frequency of vibration |
| ,: | Anhedral of the receiving point and sending point |
| : | Area of the th panel |
| : | Pressure difference on the unit length of 1/4 chord of the th panel |
| : | Steady downwash factor |
| : | Differential pressure coefficient vector |
| : | Structural mass matrix |
| : | Displacement vector |
| : | Rotor blade rotational velocity |
| : | Spline matrix |
| : | Structural mass matrix in modal space |
| : | Structural stiffness matrix at trim position |
| : | Steady structural load vector |
| : | Fundamental coupled rotating lead-lag, flap, and torsion natural frequencies, respectively |
| : | Mass radius of gyration of blade cross section () |
| : | Solidity ratio |
| : | Lift curve slope |
| : | Blade section profile drag coefficient |
| : | Elastic axis position |
| : | Cross-section moments of inertia in the flap, lead-lag, and torsion directions, respectively |
| : | Free-flow velocity |
| : | Free-flow Mach number |
| : | Normal vector of the receiving point and the sending point |
| : | Freestream velocity of the th panel |
| : | Pressure difference per unit area of the th panel |
| : | Coefficient of pressure difference |
| : | Incremental oscillatory downwash factors |
| : | Aerodynamic load vector |
| : | Structural stiffness matrix |
| : | Structural load vector |
| : | Modal matrix |
| : | Displacement vector in modal space |
| : | Structural stiffness matrix in modal space |
| : | Blade steady displacement |
| : | Steady normal downwash velocity vector |
| : | Principal mass radii of gyration of blade cross-section |
| : | Polar radius of gyration of blade cross section |
| : | Lock number |
| : | Blade precone angle |
| : | Distance to the leading edge of |
| : | Destabilizing offset between aerodynamic center and elastic axis. |
| : | Area of the blade section. |
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This study was cosupported by the National Natural Science Foundation of China (No. 11472133) and a project funded by the Priority Academic Program Development of Jiangsu Higher Education Institutions.