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

Authors: Si-wen Wang; Jing-long Han; Quan-long Chen; Hai-wei Yun; Xiao-mao Chen
Journal: International Journal of Aerospace Engineering
 (2020)

---

## 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.


---

## Body

## 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 (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 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 x,y,z induced by the surface motion can be written as follows:
(1)w¯=18π∫AU∞Δp¯Kx¯,y¯,z¯,ω,M∞qddAwhere x¯=x−ξ^,y¯=y−ζ^,z¯=z−η^, Δp¯, qd, U∞, andM∞ 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:
(2)Kx¯,y¯,z¯,ω,M∞=limnc→0∂∂nce−iωx¯/V∫−∞x¯∂∂nd1Reiωλ−M∞R/M∞dλ,where
(3)λ=x−ξ^,r=y¯2+z¯2,R=λ2+β¯2r2,β¯2=1−M∞2,ω¯=ωa∞β¯2,i=−1.Partial differentiations∂/∂nc and ∂/∂nd are taking the derivative of the receiving point and the sending point. ∂/∂nc and ∂/∂nd can be expressed by the following formulations:
(4)∂∂nc=cosnc,x∂∂x+cosnc,y∂∂y+cosnc,z∂∂z,∂∂nd=cosnd,ξ^∂∂ξ^+cosnd,ζ^∂∂ζ^+cosnd,η^∂∂η^.In the case of small perturbations, the following relations are obtained:(5)cosnd,ξ^<<cosnd,ζ^,cosnd,η^.cosnc,x cannot be ignored due to blade elastic torsion and pretwist. Formulation (4) can be rewritten as follows:
(6)∂∂nc=cosnc,x∂∂x+∂∂nc0,∂∂nd=cosnd,ζ^∂∂ζ^+cosnd,η^∂∂η^,∂∂nc0=cosnc,y∂∂y+cosnc,z∂∂z.The kernel functionK can be also written as follows:
(7)Kx¯,y¯,z¯,ω,M∞=cosnc,x−iωVe−iωx¯/V∂∂ndI¯+e−iωx¯/V∂∂x∂∂ndI¯+e−iωx¯/Vlimnc→0∂∂nc0∂∂ndI¯,where
(8)I¯=∫−∞x¯1Reiω¯λ−M∞R/M∞dλ.Because(9)∂∂ζ^1Reiω¯λ−M∞R/M∞=−∂∂y1Reiω¯λ−M∞R/M∞,∂∂η^1Reiω¯λ−M∞R/M∞=−∂∂z1Reiω¯λ−M∞R/M∞,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:
(10)cosnc,y=−sinψs,cosnc,z=cosψs,cosnd,ζ^=−sinψσ,cosnd,η^=cosψσ.And then(11)∂∂nc0∂∂nd=−sinψssinψσ∂2∂y2+cosψscosψσ∂2∂z2−sinψscosψσ+cosψssinψσ∂2∂y∂z,∂∂nd=−−sinψσ∂∂y+cosψσ∂∂z.According to the following equations:(12)∂I¯∂z=∂I¯∂r∂r∂z=z¯r∂I¯∂r,∂I¯∂y=∂I¯∂r∂r∂y=y¯r∂I¯∂r,∂2I¯∂z2=1r∂I¯∂r+z¯2r∂∂r1r∂I¯∂r,∂2I¯∂y2=1r∂I¯∂r+y¯2r∂∂r1r∂I¯∂r,∂2I¯∂yz=y¯z¯r∂∂r1r∂I¯∂r,and substituting Equations (11) and (12) into Equation (7), the kernel function can be rewritten as follows:
(13)Kx¯,y¯,z¯,ω,M∞=cosθe−iωx¯/V−sinψσy¯+cosψσz¯−iωVK1r2+∂K1∂x1r2+e−iωx¯/VK1T1r2+K2T2r4,where
(14)K1=r∂I∂r,K2=r3∂∂r1r∂I∂r,T1=cosψs−ψσ,T2=z¯cosψs−y¯sinψsz¯cosψσ−y¯sinψσ,I=−∫−∞x¯1Reiω¯λ−M∞R/M∞dλ.The computational procedure ofK1 and K2 are given by Landahl [18]. ∂K1/∂xof Equation (13) can be calculated by taking the partial derivative of K1 with respect to x. The results are given as follows:
(15)K1=I1+M∞re−ik1u1R1+u121/2,K2=−3I2−ik1M∞2r2e−ik1u1R21+u121/2−M∞re−ik1u1R1+u123/21+u12β¯2r2R2+2+M∞ru1R,∂K1∂x=∂I1∂x−ik1∂u1∂xM∞re−ik1u1R1+u121/2+−∂R∂xM∞rR21+u121/2−∂u1∂xM∞ru1R1+u123/2e−ik1u1,where
(16)u1=−x¯−M∞Rβ¯2r,R=x¯2+β¯2r2,∂u1∂x=−1β¯2r+M∞x¯β¯2rR,I0=∑n=111annc0−ik1n2c02+k12e−nc0u1,∂I0∂x=∑n=111−nc0an∂u1∂xnc0−ik1n2c02+k12e−nc0u1,J0=∑n=111ane−nc0u1n2c02+k122n2c02−k12+nc0u1n2c02+k12−ik12nc0+u1n2c02+k12.The constantsan and c0 of Equation (16) are given by Kalman et al. [19].Ifu1≥0, then I1, I2, and ∂I1/∂x will be written as follows:
(17)I1=1−u11+u121/2−ik1I0e−ik1u1,∂I1∂x=−ik1∂u1∂xI1+u121+u123/2−11+u121/2∂u1∂x−ik1∂I0∂xe−ik1u1,I2=132+ik1u11−u11+u121/2−u11+u123/2−ik1I0+k12J0e−ik1u1,And ifu1<0, then I1, I2, and ∂I1/∂x can be rewritten as follows:
(18)I1=2ReI10,k1−ReI1−u1,k1+iImI1−u1,k1,I2=2ReI20,k1−ReI2−u1,k1+iImI2−u1,k1,∂I1∂x=2Re∂I1∂x0,k1−Re∂I1∂x−u1,k1+iIm∂I1∂x−u1,k1.


#### 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 pointxi,yi,zi and the doublet point xj,yj,zj into Equations (1) and (13), the normal wash velocity w¯ij of the control point can be given as follows:
(19)w¯ij=Uj8πqdjΔp¯j∬SjKijdS.In Equation (19), Uj is the freestream velocity of the jth panel, Sj 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 qdj=0.5ρ∞Uj2, 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:
(20)Δp¯¯jlj=Δp¯jSj,where lj is the length of the 1/4 chord of the jth panel. Here, Sj=ljcosχjΔxj, where χj is the sweep angle of the doublet line and Δxj is the average chord. The pressure difference across surface Δp¯j can be written as Δp¯j=1/2ρ∞Uj2Δc¯pj, where Δc¯pj is the coefficient of the pressure difference. Substituting Sj and Δp¯j into Equation (20), Δp¯¯ can be rewritten as follows:
(21)Δp¯¯j=12ρ∞Uj2ΔxjcosχjΔc¯pj.Transform the area integral of Equation (19) into a line integral, and then the normal wash velocity is obtained as follows:
(22)w¯ij=Uj8πΔxjcosχj∫ljKijdlΔc¯pj=DijΔc¯pj,where
(23)Dij=Uj8πΔxjcosχj∫ljKijdl.The downwash factorDij can be obtained by the steady downwash factor D0ij and incremental oscillatory downwash factors D1ij and D2ij, where
(24)Dij=D0ij+D1ij+D2ij.The steady downwash factorD0ij 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 D1ij and D2ij are not presented, which can be referred to the literature [20]. The numerical form of Equation (1) in matrix notation becomes
(25)w=DΔcp,The differential pressure coefficient vectorΔcp can be obtained by Δcp=D−1w. For the jth panel, the pressure difference located at the 1/4 chord point at the midspan of the panel can be obtained byFaj=qdjΔSjΔcpj. Thereafter, the aerodynamic load vector is obtained in the form of a matrix by the following:
(26)Fa=qdΔSD−1w,where
(27)qd=qd1⋯0⋮⋱⋮0⋯qdn,ΔS=ΔS1⋯0⋮⋱⋮0⋯ΔSn.


### 2.2. Structure Modeling

Considering the geometric nonlinearity, the nonlinear forced undamped vibration equation of the blade can be simply obtained by(28)MssΩ,qsq¨s+KssΩ,qsqs=Fs,where qs is the displacement vector and Fs is the structural load vector. Mss and Kss are the structural mass matrix and stiffness matrix, respectively, and the functions of rotational speed Ω and qs. Equation (28) is transformed into the modal space by writing
(29)ps=ΦsTqs,where Φs is the modal matrix and ps 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:
(30)M¯sΩ,psp¨s+K¯sΩ,psps=F¯s,where M¯s=ΦsTMssΦs,K¯s=ΦsTKssΦs and F¯s=ΦsTFs.


### 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 matrixGas [22]. The transformation from the structural grid node motions qs into the aerodynamic grid node motions w can be given by the following:
(31)w=Gasqs.The transformation formula of the aerodynamicFa and the equivalent value acting on a structural grid Fs can be obtained by the following:
(32)Fs=GasTFa.Substituting Equations (31) and (32) into Equation (26), we can obtain the following:
(33)Fs=Qqs,where Q is a function of k and rotor rotational velocity Ω and
(34)Qk,Ω=GasTqdΔSD−1Gas.Combining Equations (28), (33), and (34), the aeroelastic equation is obtained by
(35)MssΩ,qsq¨s+KssΩ,qsqs=QΩ,qsq~sTransforming Eq. (35) to the modal space, the equation can be rewritten as follows:
(36)M¯sΩ,psp¨s+K¯sΩ,psps=Q¯k,Ωps,where Q¯k,Ω=ΦsTQk,ΩΦs.


### 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:
(37)K0ssΩ,qs0qs0=Fs0,where stiffness matrix K0ss is a function of the steady displacement qs0 and Ω and structural load vector Fs0 is given by the following:
(38)Fs0=GasTqdΔSD0−1w0.The velocity vectorw0 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 byqs00=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 theith iteration, where i=1,2,3,⋯,n, the displacements of the structural grid nodes are qs0i−1. The location of structural and aerodynamic grid nodes is refreshed. The normal downwash velocity vector w0 and the steady aerodynamic influence coefficient matrix D0−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 matrixGas is obtained using the Thin Plate Spline (TPS) interpolation method [23]. The steady structural load Fs0 is obtained solving Equation (38)
(4)
The structural grid node displacement vectorqs0i is calculated using the static structural analysis of ANSYS software when the structural load is applied to the upper wing surface
(5)
Ifqs0i−qs0i−1<ε, then, the computation converges, and if it does not, then, Steps (2)–(4) are repeatedAfter the calculation converges, the structure displacementqs0, the location of structural and aerodynamic grid nodes, and the structural load Fs0 in a static trim state are obtained. Given the large deflection and the geometric nonlinearity, the natural modes of vibration Φs0, the mass matrix M0ssΩ,qs0, and the stiffness matrix K0ssΩ,qs0 about the equilibrium position are obtained by the prestress modal analysis of ANSYS.


### 2.5. Flutter Solution

The flutter motionqs is assumed to be a small perturbation q~s about the equilibrium position qs0 and can be written as follows:
(39)qs=qs0+q~s.Substituting Equation (39) into Equation (35), we can linearize the flutter equation, based on the assumption of small perturbation, by the following:
(40)M0ssΩ,qs0q~¨s+K0ssΩ,qs0q~s+qs0=QΩ,qs0q~s+Fs0,Simplified by Equation (37), Equation (40) can be rewritten as follows:
(41)M0ssΩ,qs0q~¨s+K0ssΩ,qs0q~s=QΩ,qs0q~s.Transforming Equation (41) into the modal space, the equation can be rewritten as follows:
(42)M¯s0Ω,ps0p~¨s+K¯s0Ω,ps0p~s=Q¯k,Ωp~s,where M¯s0=Φs0TM0ssΦs0, K¯s0=Φs0TK0ssΦs0, Q¯k,Ω=Φs0TQk,ΩΦs0, p~s0=Φs0Tq~s0, and p~s=Φs0Tq~s.Based on the assumption of harmonic motion (p~s=p~sseiωt) transforming Equation (42) into a frequency domain, the aeroelastic equation can be rewritten as follows:
(43)−ω2M¯s0Ω,ps0+K¯s0Ω,ps0−Q¯k,Ωp~s=0.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 theith 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
(2)
The V-g method is used to calculate the blade flutter rotational velocityΩF and flutter frequency fF
(3)
IfΩi−ΩF<ε, the computation converges, and if it does not, Ωi+1=ΩF and Steps (1)–(2) are repeatedAfter the calculation converges, the flutter rotational velocityΩF and flutter frequency fF is obtained for the current blade pitch.


## 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 x,y,z induced by the surface motion can be written as follows:
(1)w¯=18π∫AU∞Δp¯Kx¯,y¯,z¯,ω,M∞qddAwhere x¯=x−ξ^,y¯=y−ζ^,z¯=z−η^, Δp¯, qd, U∞, andM∞ 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:
(2)Kx¯,y¯,z¯,ω,M∞=limnc→0∂∂nce−iωx¯/V∫−∞x¯∂∂nd1Reiωλ−M∞R/M∞dλ,where
(3)λ=x−ξ^,r=y¯2+z¯2,R=λ2+β¯2r2,β¯2=1−M∞2,ω¯=ωa∞β¯2,i=−1.Partial differentiations∂/∂nc and ∂/∂nd are taking the derivative of the receiving point and the sending point. ∂/∂nc and ∂/∂nd can be expressed by the following formulations:
(4)∂∂nc=cosnc,x∂∂x+cosnc,y∂∂y+cosnc,z∂∂z,∂∂nd=cosnd,ξ^∂∂ξ^+cosnd,ζ^∂∂ζ^+cosnd,η^∂∂η^.In the case of small perturbations, the following relations are obtained:(5)cosnd,ξ^<<cosnd,ζ^,cosnd,η^.cosnc,x cannot be ignored due to blade elastic torsion and pretwist. Formulation (4) can be rewritten as follows:
(6)∂∂nc=cosnc,x∂∂x+∂∂nc0,∂∂nd=cosnd,ζ^∂∂ζ^+cosnd,η^∂∂η^,∂∂nc0=cosnc,y∂∂y+cosnc,z∂∂z.The kernel functionK can be also written as follows:
(7)Kx¯,y¯,z¯,ω,M∞=cosnc,x−iωVe−iωx¯/V∂∂ndI¯+e−iωx¯/V∂∂x∂∂ndI¯+e−iωx¯/Vlimnc→0∂∂nc0∂∂ndI¯,where
(8)I¯=∫−∞x¯1Reiω¯λ−M∞R/M∞dλ.Because(9)∂∂ζ^1Reiω¯λ−M∞R/M∞=−∂∂y1Reiω¯λ−M∞R/M∞,∂∂η^1Reiω¯λ−M∞R/M∞=−∂∂z1Reiω¯λ−M∞R/M∞,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:
(10)cosnc,y=−sinψs,cosnc,z=cosψs,cosnd,ζ^=−sinψσ,cosnd,η^=cosψσ.And then(11)∂∂nc0∂∂nd=−sinψssinψσ∂2∂y2+cosψscosψσ∂2∂z2−sinψscosψσ+cosψssinψσ∂2∂y∂z,∂∂nd=−−sinψσ∂∂y+cosψσ∂∂z.According to the following equations:(12)∂I¯∂z=∂I¯∂r∂r∂z=z¯r∂I¯∂r,∂I¯∂y=∂I¯∂r∂r∂y=y¯r∂I¯∂r,∂2I¯∂z2=1r∂I¯∂r+z¯2r∂∂r1r∂I¯∂r,∂2I¯∂y2=1r∂I¯∂r+y¯2r∂∂r1r∂I¯∂r,∂2I¯∂yz=y¯z¯r∂∂r1r∂I¯∂r,and substituting Equations (11) and (12) into Equation (7), the kernel function can be rewritten as follows:
(13)Kx¯,y¯,z¯,ω,M∞=cosθe−iωx¯/V−sinψσy¯+cosψσz¯−iωVK1r2+∂K1∂x1r2+e−iωx¯/VK1T1r2+K2T2r4,where
(14)K1=r∂I∂r,K2=r3∂∂r1r∂I∂r,T1=cosψs−ψσ,T2=z¯cosψs−y¯sinψsz¯cosψσ−y¯sinψσ,I=−∫−∞x¯1Reiω¯λ−M∞R/M∞dλ.The computational procedure ofK1 and K2 are given by Landahl [18]. ∂K1/∂xof Equation (13) can be calculated by taking the partial derivative of K1 with respect to x. The results are given as follows:
(15)K1=I1+M∞re−ik1u1R1+u121/2,K2=−3I2−ik1M∞2r2e−ik1u1R21+u121/2−M∞re−ik1u1R1+u123/21+u12β¯2r2R2+2+M∞ru1R,∂K1∂x=∂I1∂x−ik1∂u1∂xM∞re−ik1u1R1+u121/2+−∂R∂xM∞rR21+u121/2−∂u1∂xM∞ru1R1+u123/2e−ik1u1,where
(16)u1=−x¯−M∞Rβ¯2r,R=x¯2+β¯2r2,∂u1∂x=−1β¯2r+M∞x¯β¯2rR,I0=∑n=111annc0−ik1n2c02+k12e−nc0u1,∂I0∂x=∑n=111−nc0an∂u1∂xnc0−ik1n2c02+k12e−nc0u1,J0=∑n=111ane−nc0u1n2c02+k122n2c02−k12+nc0u1n2c02+k12−ik12nc0+u1n2c02+k12.The constantsan and c0 of Equation (16) are given by Kalman et al. [19].Ifu1≥0, then I1, I2, and ∂I1/∂x will be written as follows:
(17)I1=1−u11+u121/2−ik1I0e−ik1u1,∂I1∂x=−ik1∂u1∂xI1+u121+u123/2−11+u121/2∂u1∂x−ik1∂I0∂xe−ik1u1,I2=132+ik1u11−u11+u121/2−u11+u123/2−ik1I0+k12J0e−ik1u1,And ifu1<0, then I1, I2, and ∂I1/∂x can be rewritten as follows:
(18)I1=2ReI10,k1−ReI1−u1,k1+iImI1−u1,k1,I2=2ReI20,k1−ReI2−u1,k1+iImI2−u1,k1,∂I1∂x=2Re∂I1∂x0,k1−Re∂I1∂x−u1,k1+iIm∂I1∂x−u1,k1.


### 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 pointxi,yi,zi and the doublet point xj,yj,zj into Equations (1) and (13), the normal wash velocity w¯ij of the control point can be given as follows:
(19)w¯ij=Uj8πqdjΔp¯j∬SjKijdS.In Equation (19), Uj is the freestream velocity of the jth panel, Sj 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 qdj=0.5ρ∞Uj2, 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:
(20)Δp¯¯jlj=Δp¯jSj,where lj is the length of the 1/4 chord of the jth panel. Here, Sj=ljcosχjΔxj, where χj is the sweep angle of the doublet line and Δxj is the average chord. The pressure difference across surface Δp¯j can be written as Δp¯j=1/2ρ∞Uj2Δc¯pj, where Δc¯pj is the coefficient of the pressure difference. Substituting Sj and Δp¯j into Equation (20), Δp¯¯ can be rewritten as follows:
(21)Δp¯¯j=12ρ∞Uj2ΔxjcosχjΔc¯pj.Transform the area integral of Equation (19) into a line integral, and then the normal wash velocity is obtained as follows:
(22)w¯ij=Uj8πΔxjcosχj∫ljKijdlΔc¯pj=DijΔc¯pj,where
(23)Dij=Uj8πΔxjcosχj∫ljKijdl.The downwash factorDij can be obtained by the steady downwash factor D0ij and incremental oscillatory downwash factors D1ij and D2ij, where
(24)Dij=D0ij+D1ij+D2ij.The steady downwash factorD0ij 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 D1ij and D2ij are not presented, which can be referred to the literature [20]. The numerical form of Equation (1) in matrix notation becomes
(25)w=DΔcp,The differential pressure coefficient vectorΔcp can be obtained by Δcp=D−1w. For the jth panel, the pressure difference located at the 1/4 chord point at the midspan of the panel can be obtained byFaj=qdjΔSjΔcpj. Thereafter, the aerodynamic load vector is obtained in the form of a matrix by the following:
(26)Fa=qdΔSD−1w,where
(27)qd=qd1⋯0⋮⋱⋮0⋯qdn,ΔS=ΔS1⋯0⋮⋱⋮0⋯ΔSn.


## 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 x,y,z induced by the surface motion can be written as follows:
(1)w¯=18π∫AU∞Δp¯Kx¯,y¯,z¯,ω,M∞qddAwhere x¯=x−ξ^,y¯=y−ζ^,z¯=z−η^, Δp¯, qd, U∞, andM∞ 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:
(2)Kx¯,y¯,z¯,ω,M∞=limnc→0∂∂nce−iωx¯/V∫−∞x¯∂∂nd1Reiωλ−M∞R/M∞dλ,where
(3)λ=x−ξ^,r=y¯2+z¯2,R=λ2+β¯2r2,β¯2=1−M∞2,ω¯=ωa∞β¯2,i=−1.Partial differentiations∂/∂nc and ∂/∂nd are taking the derivative of the receiving point and the sending point. ∂/∂nc and ∂/∂nd can be expressed by the following formulations:
(4)∂∂nc=cosnc,x∂∂x+cosnc,y∂∂y+cosnc,z∂∂z,∂∂nd=cosnd,ξ^∂∂ξ^+cosnd,ζ^∂∂ζ^+cosnd,η^∂∂η^.In the case of small perturbations, the following relations are obtained:(5)cosnd,ξ^<<cosnd,ζ^,cosnd,η^.cosnc,x cannot be ignored due to blade elastic torsion and pretwist. Formulation (4) can be rewritten as follows:
(6)∂∂nc=cosnc,x∂∂x+∂∂nc0,∂∂nd=cosnd,ζ^∂∂ζ^+cosnd,η^∂∂η^,∂∂nc0=cosnc,y∂∂y+cosnc,z∂∂z.The kernel functionK can be also written as follows:
(7)Kx¯,y¯,z¯,ω,M∞=cosnc,x−iωVe−iωx¯/V∂∂ndI¯+e−iωx¯/V∂∂x∂∂ndI¯+e−iωx¯/Vlimnc→0∂∂nc0∂∂ndI¯,where
(8)I¯=∫−∞x¯1Reiω¯λ−M∞R/M∞dλ.Because(9)∂∂ζ^1Reiω¯λ−M∞R/M∞=−∂∂y1Reiω¯λ−M∞R/M∞,∂∂η^1Reiω¯λ−M∞R/M∞=−∂∂z1Reiω¯λ−M∞R/M∞,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:
(10)cosnc,y=−sinψs,cosnc,z=cosψs,cosnd,ζ^=−sinψσ,cosnd,η^=cosψσ.And then(11)∂∂nc0∂∂nd=−sinψssinψσ∂2∂y2+cosψscosψσ∂2∂z2−sinψscosψσ+cosψssinψσ∂2∂y∂z,∂∂nd=−−sinψσ∂∂y+cosψσ∂∂z.According to the following equations:(12)∂I¯∂z=∂I¯∂r∂r∂z=z¯r∂I¯∂r,∂I¯∂y=∂I¯∂r∂r∂y=y¯r∂I¯∂r,∂2I¯∂z2=1r∂I¯∂r+z¯2r∂∂r1r∂I¯∂r,∂2I¯∂y2=1r∂I¯∂r+y¯2r∂∂r1r∂I¯∂r,∂2I¯∂yz=y¯z¯r∂∂r1r∂I¯∂r,and substituting Equations (11) and (12) into Equation (7), the kernel function can be rewritten as follows:
(13)Kx¯,y¯,z¯,ω,M∞=cosθe−iωx¯/V−sinψσy¯+cosψσz¯−iωVK1r2+∂K1∂x1r2+e−iωx¯/VK1T1r2+K2T2r4,where
(14)K1=r∂I∂r,K2=r3∂∂r1r∂I∂r,T1=cosψs−ψσ,T2=z¯cosψs−y¯sinψsz¯cosψσ−y¯sinψσ,I=−∫−∞x¯1Reiω¯λ−M∞R/M∞dλ.The computational procedure ofK1 and K2 are given by Landahl [18]. ∂K1/∂xof Equation (13) can be calculated by taking the partial derivative of K1 with respect to x. The results are given as follows:
(15)K1=I1+M∞re−ik1u1R1+u121/2,K2=−3I2−ik1M∞2r2e−ik1u1R21+u121/2−M∞re−ik1u1R1+u123/21+u12β¯2r2R2+2+M∞ru1R,∂K1∂x=∂I1∂x−ik1∂u1∂xM∞re−ik1u1R1+u121/2+−∂R∂xM∞rR21+u121/2−∂u1∂xM∞ru1R1+u123/2e−ik1u1,where
(16)u1=−x¯−M∞Rβ¯2r,R=x¯2+β¯2r2,∂u1∂x=−1β¯2r+M∞x¯β¯2rR,I0=∑n=111annc0−ik1n2c02+k12e−nc0u1,∂I0∂x=∑n=111−nc0an∂u1∂xnc0−ik1n2c02+k12e−nc0u1,J0=∑n=111ane−nc0u1n2c02+k122n2c02−k12+nc0u1n2c02+k12−ik12nc0+u1n2c02+k12.The constantsan and c0 of Equation (16) are given by Kalman et al. [19].Ifu1≥0, then I1, I2, and ∂I1/∂x will be written as follows:
(17)I1=1−u11+u121/2−ik1I0e−ik1u1,∂I1∂x=−ik1∂u1∂xI1+u121+u123/2−11+u121/2∂u1∂x−ik1∂I0∂xe−ik1u1,I2=132+ik1u11−u11+u121/2−u11+u123/2−ik1I0+k12J0e−ik1u1,And ifu1<0, then I1, I2, and ∂I1/∂x can be rewritten as follows:
(18)I1=2ReI10,k1−ReI1−u1,k1+iImI1−u1,k1,I2=2ReI20,k1−ReI2−u1,k1+iImI2−u1,k1,∂I1∂x=2Re∂I1∂x0,k1−Re∂I1∂x−u1,k1+iIm∂I1∂x−u1,k1.


## 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 pointxi,yi,zi and the doublet point xj,yj,zj into Equations (1) and (13), the normal wash velocity w¯ij of the control point can be given as follows:
(19)w¯ij=Uj8πqdjΔp¯j∬SjKijdS.In Equation (19), Uj is the freestream velocity of the jth panel, Sj 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 qdj=0.5ρ∞Uj2, 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:
(20)Δp¯¯jlj=Δp¯jSj,where lj is the length of the 1/4 chord of the jth panel. Here, Sj=ljcosχjΔxj, where χj is the sweep angle of the doublet line and Δxj is the average chord. The pressure difference across surface Δp¯j can be written as Δp¯j=1/2ρ∞Uj2Δc¯pj, where Δc¯pj is the coefficient of the pressure difference. Substituting Sj and Δp¯j into Equation (20), Δp¯¯ can be rewritten as follows:
(21)Δp¯¯j=12ρ∞Uj2ΔxjcosχjΔc¯pj.Transform the area integral of Equation (19) into a line integral, and then the normal wash velocity is obtained as follows:
(22)w¯ij=Uj8πΔxjcosχj∫ljKijdlΔc¯pj=DijΔc¯pj,where
(23)Dij=Uj8πΔxjcosχj∫ljKijdl.The downwash factorDij can be obtained by the steady downwash factor D0ij and incremental oscillatory downwash factors D1ij and D2ij, where
(24)Dij=D0ij+D1ij+D2ij.The steady downwash factorD0ij 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 D1ij and D2ij are not presented, which can be referred to the literature [20]. The numerical form of Equation (1) in matrix notation becomes
(25)w=DΔcp,The differential pressure coefficient vectorΔcp can be obtained by Δcp=D−1w. For the jth panel, the pressure difference located at the 1/4 chord point at the midspan of the panel can be obtained byFaj=qdjΔSjΔcpj. Thereafter, the aerodynamic load vector is obtained in the form of a matrix by the following:
(26)Fa=qdΔSD−1w,where
(27)qd=qd1⋯0⋮⋱⋮0⋯qdn,ΔS=ΔS1⋯0⋮⋱⋮0⋯ΔSn.


## 2.2. Structure Modeling

Considering the geometric nonlinearity, the nonlinear forced undamped vibration equation of the blade can be simply obtained by(28)MssΩ,qsq¨s+KssΩ,qsqs=Fs,where qs is the displacement vector and Fs is the structural load vector. Mss and Kss are the structural mass matrix and stiffness matrix, respectively, and the functions of rotational speed Ω and qs. Equation (28) is transformed into the modal space by writing
(29)ps=ΦsTqs,where Φs is the modal matrix and ps 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:
(30)M¯sΩ,psp¨s+K¯sΩ,psps=F¯s,where M¯s=ΦsTMssΦs,K¯s=ΦsTKssΦs and F¯s=ΦsTFs.


## 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 matrixGas [22]. The transformation from the structural grid node motions qs into the aerodynamic grid node motions w can be given by the following:
(31)w=Gasqs.The transformation formula of the aerodynamicFa and the equivalent value acting on a structural grid Fs can be obtained by the following:
(32)Fs=GasTFa.Substituting Equations (31) and (32) into Equation (26), we can obtain the following:
(33)Fs=Qqs,where Q is a function of k and rotor rotational velocity Ω and
(34)Qk,Ω=GasTqdΔSD−1Gas.Combining Equations (28), (33), and (34), the aeroelastic equation is obtained by
(35)MssΩ,qsq¨s+KssΩ,qsqs=QΩ,qsq~sTransforming Eq. (35) to the modal space, the equation can be rewritten as follows:
(36)M¯sΩ,psp¨s+K¯sΩ,psps=Q¯k,Ωps,where Q¯k,Ω=ΦsTQk,ΩΦs.


## 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:
(37)K0ssΩ,qs0qs0=Fs0,where stiffness matrix K0ss is a function of the steady displacement qs0 and Ω and structural load vector Fs0 is given by the following:
(38)Fs0=GasTqdΔSD0−1w0.The velocity vectorw0 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 byqs00=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 theith iteration, where i=1,2,3,⋯,n, the displacements of the structural grid nodes are qs0i−1. The location of structural and aerodynamic grid nodes is refreshed. The normal downwash velocity vector w0 and the steady aerodynamic influence coefficient matrix D0−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 matrixGas is obtained using the Thin Plate Spline (TPS) interpolation method [23]. The steady structural load Fs0 is obtained solving Equation (38)
(4)
The structural grid node displacement vectorqs0i is calculated using the static structural analysis of ANSYS software when the structural load is applied to the upper wing surface
(5)
Ifqs0i−qs0i−1<ε, then, the computation converges, and if it does not, then, Steps (2)–(4) are repeatedAfter the calculation converges, the structure displacementqs0, the location of structural and aerodynamic grid nodes, and the structural load Fs0 in a static trim state are obtained. Given the large deflection and the geometric nonlinearity, the natural modes of vibration Φs0, the mass matrix M0ssΩ,qs0, and the stiffness matrix K0ssΩ,qs0 about the equilibrium position are obtained by the prestress modal analysis of ANSYS.


## 2.5. Flutter Solution

The flutter motionqs is assumed to be a small perturbation q~s about the equilibrium position qs0 and can be written as follows:
(39)qs=qs0+q~s.Substituting Equation (39) into Equation (35), we can linearize the flutter equation, based on the assumption of small perturbation, by the following:
(40)M0ssΩ,qs0q~¨s+K0ssΩ,qs0q~s+qs0=QΩ,qs0q~s+Fs0,Simplified by Equation (37), Equation (40) can be rewritten as follows:
(41)M0ssΩ,qs0q~¨s+K0ssΩ,qs0q~s=QΩ,qs0q~s.Transforming Equation (41) into the modal space, the equation can be rewritten as follows:
(42)M¯s0Ω,ps0p~¨s+K¯s0Ω,ps0p~s=Q¯k,Ωp~s,where M¯s0=Φs0TM0ssΦs0, K¯s0=Φs0TK0ssΦs0, Q¯k,Ω=Φs0TQk,ΩΦs0, p~s0=Φs0Tq~s0, and p~s=Φs0Tq~s.Based on the assumption of harmonic motion (p~s=p~sseiωt) transforming Equation (42) into a frequency domain, the aeroelastic equation can be rewritten as follows:
(43)−ω2M¯s0Ω,ps0+K¯s0Ω,ps0−Q¯k,Ωp~s=0.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 theith 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
(2)
The V-g method is used to calculate the blade flutter rotational velocityΩF and flutter frequency fF
(3)
IfΩi−ΩF<ε, the computation converges, and if it does not, Ωi+1=ΩF and Steps (1)–(2) are repeatedAfter the calculation converges, the flutter rotational velocityΩF and flutter frequency fF 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ωw=1.15, ωv=1.5, and ωϕ=3.0 or 5.0, respectively. Table 1 shows other blade nondimensional parameters [3].Table 1
Values of parameters of uniform blade for numerical results.
km1/R=0
c/R=π/40
a=6
km2/R=0.025
σ=0.1
βp=0.05
kA/km=1.5
γ=5
CD0=0.0095Figure1 exhibits the steady tip deflection wtip, flap and lead-lag bending deflection vtip, 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 1
Comparison of steady tip deflection by the MDLM and the quasisteady strip theory (ωw=1.15, ωv=1.5, and ωϕ=5.0).
(a)
(b)
(c)Figure2 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.Figure 2
Effect of precone on the stability boundaries (ωw=1.15, ωv=1.5, and ωϕ==3.0).


### 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.Figure 3
Meshing of the blade model.Figure 4
Natural frequencies of the blade (Λ=45°).Figure 5
Natural frequencies of the blade (Λ=15°).


### 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 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.Figure 6
Blade section.Table 2
Values of parameters of blade section.
A=0.001229m2
Iz=0.441×10−5m4
XA=0.006704m
Iy=0.766×10−7m4
J=0.244×10−6m4
c=0.215m
XE=0.060454mThe isotropic material is aluminum. Figure7 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.367m 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.Figure 7
Blade planform with tip sweep (Λ=45°).Table3 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.Table 3
Natural frequency of hingeless rotor blade (AR=18.07).
Number of elements
Mode 1 (Hz)
Mode 2 (Hz)
Mode 3 (Hz)
Mode 4 (Hz)
4880
7.57
11.837
15.508
26.472
5490
7.57
11.837
15.506
26.471
6100
7.5697
11.836
15.505
26.469
6710
7.5696
11.836
15.504
26.467Figures8–9 show the flutter rotational velocities, and Figures 10-11 show the V-g diagrams, with respect to various aspect ratios, where Ω=500rpm. 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 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%.Figure 8
Flutter rotational velocities of the blade (AR=18.07).Figure 9
Flutter rotational velocities of the blade (AR=12.5).Figure 10
V-g diagram of the blade (AR=18.07).Figure 11
V-g diagram of the blade (AR=12.5).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 (theta0=6°). Figures 14-15 show the shapes of these modes of the blade (AR=12.5). 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. Figures16–17 show the displacements of the 3-D finite element nodes at r=0.75R 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 ratioFigure 12
Natural frequencies (AR=18.07).Figure 13
Natural frequencies (AR=12.5).Figure 14
Lead-lag/torsion coupled mode.
(a)
Front view
(b)
Top viewFigure 15
Flap/torsion coupled mode.Figure 16
Deformed sections of 3-D model (AR=18.07).Figure 17
Deformed section of 3-D model (AR=12.5).


## 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ωw=1.15, ωv=1.5, and ωϕ=3.0 or 5.0, respectively. Table 1 shows other blade nondimensional parameters [3].Table 1
Values of parameters of uniform blade for numerical results.
km1/R=0
c/R=π/40
a=6
km2/R=0.025
σ=0.1
βp=0.05
kA/km=1.5
γ=5
CD0=0.0095Figure1 exhibits the steady tip deflection wtip, flap and lead-lag bending deflection vtip, 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 1
Comparison of steady tip deflection by the MDLM and the quasisteady strip theory (ωw=1.15, ωv=1.5, and ωϕ=5.0).
(a)
(b)
(c)Figure2 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.Figure 2
Effect of precone on the stability boundaries (ωw=1.15, ωv=1.5, and ωϕ==3.0).


## 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.Figure 3
Meshing of the blade model.Figure 4
Natural frequencies of the blade (Λ=45°).Figure 5
Natural frequencies of the blade (Λ=15°).


## 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 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.Figure 6
Blade section.Table 2
Values of parameters of blade section.
A=0.001229m2
Iz=0.441×10−5m4
XA=0.006704m
Iy=0.766×10−7m4
J=0.244×10−6m4
c=0.215m
XE=0.060454mThe isotropic material is aluminum. Figure7 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.367m 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.Figure 7
Blade planform with tip sweep (Λ=45°).Table3 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.Table 3
Natural frequency of hingeless rotor blade (AR=18.07).
Number of elements
Mode 1 (Hz)
Mode 2 (Hz)
Mode 3 (Hz)
Mode 4 (Hz)
4880
7.57
11.837
15.508
26.472
5490
7.57
11.837
15.506
26.471
6100
7.5697
11.836
15.505
26.469
6710
7.5696
11.836
15.504
26.467Figures8–9 show the flutter rotational velocities, and Figures 10-11 show the V-g diagrams, with respect to various aspect ratios, where Ω=500rpm. 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 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%.Figure 8
Flutter rotational velocities of the blade (AR=18.07).Figure 9
Flutter rotational velocities of the blade (AR=12.5).Figure 10
V-g diagram of the blade (AR=18.07).Figure 11
V-g diagram of the blade (AR=12.5).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 (theta0=6°). Figures 14-15 show the shapes of these modes of the blade (AR=12.5). 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. Figures16–17 show the displacements of the 3-D finite element nodes at r=0.75R 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 ratioFigure 12
Natural frequencies (AR=18.07).Figure 13
Natural frequencies (AR=12.5).Figure 14
Lead-lag/torsion coupled mode.
(a)
Front view
(b)
Top viewFigure 15
Flap/torsion coupled mode.Figure 16
Deformed sections of 3-D model (AR=18.07).Figure 17
Deformed section of 3-D model (AR=12.5).


## 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



---

*Source: 1891765-2020-02-17.xml*