Free Vibration Analysis of Rotating Pretwisted Functionally Graded Sandwich Blades

A new structural dynamic model for the free vibration characteristic analysis of rotating pretwisted functionally graded (FG) sandwich blades is developed. The sandwich blade is made up of two functionally graded skins and a homogeneous material core. The thick shell theory is applied to derive the basic equations of motion of the rotating FG sandwich blade by considering the effects of centrifugal and Coriolis forces. The mode shapes are expanded in terms of two-dimensional algebraic polynomials in the Rayleigh–Ritz method, and the static and dynamic natural frequencies of the blade are obtained. The convergence analysis is studied, and the accuracy of the proposed model is verified by comparing with the literature results and ANSYS data. The effects of frequency parameters such as the twist angle, the thickness ratio, the aspect ratio, the layer thickness ratio, the scalar parameter of volume fraction, the stagger angle, and the rotation velocity on the vibration characteristics for pretwist FG sandwich blade are investigated in detail. In addition, the phenomena of frequency locus veering and mode shape exchanging occur in the static and dynamic states. Frequency locus veering is essentially caused by the coupling between different modes.


Introduction
The aero gas turbine engine is the key technology of the aircraft; blade is one of the most important parts of the system.Vibration, especially resonance of the blade, will produce larger stress, which will lead to fatigue failure.Thus, the dynamic behavior of rotating blades has been studied by numerous researchers.
Some simplified blade models such as beams, plates, and shells are used to investigate the vibration characteristics of the blades in the published literatures.Leissa and Jacob [1] investigated the free vibration of pretwisted, cantilevered beams and plates by using the Ritz method.The work was the first three-dimensional study of the problem.Yoo et al. [2] derived the equations of motion with a concentrated mass by using a modeling method with hybrid deformation variables and investigated the effects of the nondimensional parameters on the vibration characteristics of the pretwisted rotating blade through numerical analysis.Chandiramani et al. [3] simplified the pretwisted rotating blade as a laminated composite, hollow uniform box-beam model, which considers the centrifugal and Coriolis effects, transverse shear flexibility, and restrained warping, and studied the free and forced vibration by using HSDT and Galerkin method.Carrera et al. [4] used the Carrera unified formulation and FEM to study the free vibration analysis of rotating blades; the Coriolis and centrifugal force fields were included in their work.
The beam models are quite suitable for blades with large aspect ratio and low width thickness ratio.However, flexible rotating structures such as blades with low aspect ratios are also widely used in actual engineering applications.Therefore, more and more literatures were found to study the vibration of blade with plate and shell models.
Qatu and Leissa [5] were the first to study the effect of plate parameters such as twist angle, stacking sequence, and lamination angle on the natural frequency and mode shapes of laminated composite pretwisted cantilever plates by using the shallow shell theory and Ritz method.Nabi and Ganesan [6] analyzed the vibration characteristics of metal matrix composite pretwisted blades by using beam and plate theories, respectively, and summed up the quantitative comparison of natural frequency.Yoo and Chung [7] developed a linear dynamic modeling method for plates by using two in-plane stretch variables and analyzed and studied the transient characteristics of rotating plates.Hu et al. [8] proposed a numerical procedure for the free vibration analysis of pretwisted thin plates based on the thin shell theory and studied the vibration characteristics considering different twist rates and aspect ratios.Hashemi et al. [9] used the Mindlin plate theory and Kane dynamic method to develop a finite element formulation for vibration analysis of rotating thick plates.The coupling between in-plane and out-of-plane deformations and Coriolis effect was considered in the study.Sinha and Turner [10] derived the governing partial differential equation of motion for the rotating pretwisted plate by using the thin shell theory and studied the free vibration of a turbomachinery cantilevered airfoil blade with the Rayleigh-Ritz technique by considering the centrifugal force filed as a quasi-static load.Sun et al. [11] presented a dynamic model based on CLPT, investigated the vibration behavior of a rotating blade by using Hamilton's principle, and studied the point and distribution forced response using a proportional damping model.
Rao and Gupta [12] studied the vibration of natural frequencies with various parameters of a rotating pretwisted blade with a small aspect ratio by using the classical bending theory of thin shells.Sun et al. [13] developed a novel dynamic model for pretwisted rotating compressor blades by using the general shell theory and studied the eigenfrequencies and damping properties of the pretwisted rotating blade.Sinha and Zylka [14] derived the governing partial differential equation of motion based on the thin shell theory and formulated the free vibration of a turbomachinery cantilevered airfoil by considering it as an anisotropic shell in a centrifugal force field.Kielb et al. [15] gave the complete theoretical and experimental results of a joint research study on the vibration characteristics of pretwisted cantilever plates.
Claassen and Thorne [16] were the first to find out and discuss the curve veering in the plate vibration based on CLPT.Afolabi and Mehmed [17] studied the curve veering and flutter of rotating blades based on the theory of algebraic curves and catastrophe theory.They found that the frequency loci of undamped rotating blades do not cross, but must either repel each other or attract each other.Yoo and Pierre [18,19] studied the vibration characteristics of rotating cantilever rectangular plates and composite plates, respectively, and found out and discussed the natural frequency locus veering, crossing, and associated mode shape variations in detail.
Most of the early studies were focused on homogeneous materials, and in recent years, some researchers have extended the range to composite materials because of their excellent properties.The structure was developed from single layer to multilayer and composite layer, and more and more literature has been published.
Oh et al. [20] studied the vibration of turbomachinery rotating blades made up of functionally graded materials and investigated the influence of parameters on the frequency of rotating blades.A refined dynamic theory involving the coupling between flapping, lagging, and transverse shear was used in their work.Ramesh and Rao [21] studied the natural frequencies of a pretwisted rotating FG cantilever beam by using Lagrange's equation and the Rayleigh-Ritz method.The influences of different parameters and coupling between chordwise and flapwise bending modes on the natural frequency were investigated.Oh and Yoo [22] presented a new structural dynamic model to study the vibration characteristics of rotating blades.The blade was modeled as a functionally graded material pretwisted beam by using the Rayleigh-Ritz and Kane method.Mantari and Granados [23] studied the free vibration of FG plates by using a new novel FSDT with 4 unknowns.Li and Zhang [24] developed a dynamic model for FG plates undergoing large overall motions, studied the free vibration of rotating cantilever FGM rectangular plates, and found frequency locus veering and associated mode shift phenomena.Sun et al. [25] developed a novel dynamic model for multilayer blades by using the quadratic layerwise theory and studied the structural dynamics, particularly the damping properties of the rotating blade.Frequency locus veering with rotation speed was also found.Cao et al. [26] developed a rotating cantilever pretwisted sandwich plate model with a prestagger angle and investigated the vibrational behavior of a turbine blade with thermal barrier coating layers by using FSDT, von Karman plate theory, and Chebyshev-Ritz method.
In the present paper, a new structural dynamic model is developed to study the free vibration characteristics of pretwisted rotating FG sandwich blades.The thick shell theory is applied to derive the basic equations of motion by considering the effects of centrifugal and Coriolis forces.The accuracy of the proposed model is verified by comparing with the literature and ANSYS.The effects of frequency parameters such as the twist angle, the thickness ratio, the aspect ratio, the layer thickness ratio, the scalar parameter of volume fraction, the stagger angle, and the rotation velocity on the vibration characteristics are investigated in detail.Furthermore, frequency locus veering and mode shape exchanging phenomena are found both in both static and dynamic states.

Theoretical Analysis
2.1.Basic Equations.As shown in Figure 1, the functionally graded sandwich blade is modeled as a cantilever pretwisted thick plate with the twist angle θ at the free end, which is clamped to the rigid disk with a radius R, mounted with a stagger angle φ.The geometric parameters of the blade are the span (length) dimension L, the chord (width) dimension b, the total thickness h, and the angular rotating velocity about the rigid disk axis Ω.
In this paper, three coordinate systems are established for dynamic modeling.One is the XYZ-coordinate system (with unit vectors (i X , i Y , i Z )), where the x-axis is along the spanwise direction of the blade, Y is the rotation axis, and the Z-axis is perpendicular to the XY-plane following the right-hand rule.The other is the xyz-coordinate system (with unit vectors (i x , i y , i z )), where the origin lies in the root of the mid-surface of the blade and the x-axis is parallel to the X-axis; the y-axis can be obtained by rotating around the x-axis with the stagger angle φ starting from the Y-axis International Journal of Aerospace Engineering direction.Figure 2 shows the last xη coordinate, a dynamic surface coordinate where the η-axis is perpendicular to the x-axis and equal to the y-axis at x = 0, rotating at a uniform twist rate k = θ/L.Thus, the rotating angle at x can be expressed as The position vector r 0 0 x, η of a typical point on the mid-surface of the pretwisted blade in the xη-coordinate can be written as Based on the Frenet-Serret formula [28], the corresponding unit base vectors can be expressed as where A is the Lame parameter of the middle surface in the x-direction, given by Then, an arbitrary point of the blade before deformation can be expressed by a position vector r 0 as As illustrated in Figure 3, the core of the blade (h 1 ≤ ζ ≤ h 2 ) is fully metal (isotropic) and two skins are composed of functionally graded material (FGM) in the thickness direction.The top skin varies from ceramic-rich surface (ζ = h 3 ) to metal-rich surface (ζ = h 2 ), while the bottom skin varies from ceramic-rich surface (ζ = h 0 ) to metal-rich surface (ζ = h 1 ).There are no interfaces between the core and skins of the functionally graded sandwich blade, and two skins are symmetrical about the middle As shown in Figure 4, n is the scalar parameter that allows users to define the gradation of material properties in the thickness direction.The n = 0 case corresponds to a fully metal blade.The volume fraction for the ceramic phase is given as In the present work, the homogenization procedure and the law of mixtures are used.The variation of the parameters in the thickness direction is given by where P ζ is the physical parameter, which refers to density ρ, elastic modulus E, and Poisson's ratio μ.In this paper, P m and P c are the physical parameters of metal and ceramic phase.

The Strain Energy.
According to the first-order shear deformation theory, the displacements can be written as where ϕ x and ϕ η are mid-surface rotations and u 0 , v 0 , and w 0 are mid-surface displacements of the blade along the x, η, and ζ directions.
According to the thick shell theory [29] and the improved version of the Novozhilov nonlinear shell theory [30], the strains at any point in the blade can be written as  4 International Journal of Aerospace Engineering where where R xη is the radius of the twist, which can be derived as [31] 1 The nonlinear parts, which are underlined in (11), will be ignored in the free vibration analysis of the blade.
The stress-strain equations can be written as where the elastic constants Q ij ζ are functions of thickness ζ, which is defined as where K s is the shear correction coefficient.Since the shear correction coefficient depends on the blade material properties and boundary conditions, it is difficult to obtain the exact value.In this paper, the shear correction coefficient K s will be selected uniformly as 5/6 [32].
The strain energy of the blade is defined as where A ij , B ij , and D ij are the extensional, bending-extensional coupling, and bending stiffness [32], A ij ′ and B ij ′ are the stiffness due to the initial twist of the blade.They are defined as 16b It should be noted that the stiffness B ij and B ij ′ will be zero while the structure is symmetrical about the middle surface in the thickness direction.Figure 5 shows the variations of the Then, the angular velocity ω in the xyz-coordinate of the pretwisted plate is expressed as The centrifugal force F C per unit volume in the blade is given by where Φ is the sum of the stagger angle and the twist angle of the point, given by The rotation of the blade about the turbine axial direction produces a significant amount of membrane stresses, which are essentially acting in the radial direction (the spanwise and chordwise directions of the blade).The component of centrifugal force along the spanwise direction σ cs and the chordwise direction σ cc can be written as [10,14] where Θ is the sum of the stagger angle and the total twist angle, given by The centrifugal potential energy can be derived as According to [26], ε cs and ε cc are given by Substituting ( 21) and ( 24) into ( 23), the potential energy can be written as

25
The total potential energy can be written as  An arbitrary point of the blade after deformation can be expressed by the position vector r as The velocity V r of the point due to the rotation can be derived as The kinetic energy reads where where "•" is the first-order derivation of time; I 0 , I 1 , and I 2 in ( 25) and (30) are defined in terms of the density ρ ζ as

Frequency Solving
The Rayleigh-Ritz method is used to solve the natural frequency of the blade; the expressions of the displacements and rotations are assumed The functions U N x, η , V N x, η , W N x, η , Φ xN x, η , and Φ ηN x, η represent the N order mode shape, and ω N is the corresponding circular frequency.The algebraic polynomials can form a mathematically complete set of functions, which guarantees convergence to the exact solution as the number of terms taken increases [5].Thus, in this paper, the mode shapes are expanded in terms of two-dimensional algebraic polynomials in terms of the nondimensional coordinates χ and ξ as

9
International Journal of Aerospace Engineering W mn χ m ξ n , 33c where χ = x/L and ξ = 2η/b.U ij , V kl , W mn , Φ xij , and Φ ηkl are undetermined coefficients.The Rayleigh-Ritz method imposes a necessary condition that trial functions must satisfy the geometric boundary conditions arbitrarily.The cantilever blade is clamped at χ = 0 (u = v = w = dw/dx = 0) along the spanwise direction and is completely free at χ = 1, ξ = −1, and ξ = 1.The indices i, k, and m in (33), (46), and (47) begin with i = 1, k = 1, and m = 2; it can be ensured that all terms of polynomials satisfied the geometric boundary conditions.
To solve the free frequency of the blade, ( 32) and ( 33) are substituted into (26) and (30) in order to obtain the maximum strain energy (U max ) and kinetic energy (T max ) by setting t = 0.The Rayleigh-Ritz method requires the minimization of the function = T max − U max , which can be achieved by taking the derivatives There are a total of 2 × I × J + 1 + 2 × K × L + 1 + M − 1 × N + 1 equations, which can be described by where K and M are the stiffness and mass matrices, respectively; Δ is the vector of undetermined coefficients, given by The generalized eigenvalues (natural frequency) can be obtained by setting the coefficient matrix of (35) equal to zero, and the corresponding mode shape can be determined by substituting the eigenvector Δ back into (33).m /E m for the functionally graded sandwich blade in the following section, and the nondimensional rotation velocity is Ω = Ω/ω 1 , where ω 1 is the first static circular frequency (Ω = 0) of the functionally graded sandwich blade.
According to [5], the number of terms chosen for each component of the displacements is the same.The results of the homogeneous material blade and functionally graded sandwich blade are given in Tables 1 and 2. The maximum difference between the 5 × 9 × 9 = 405 term and the 5 × 8 × 8 = 320 term solutions for the first six modes of the homogeneous material blade is no more than 0.28%, and 0.21% for the FG sandwich blade.But the time spent in computing has multiplied.Moreover, the eigenvalues are more likely to be ill-conditioned with the increase in the size of the matrix in numerical calculation.Therefore, it has been decided to use 5 × 8 × 8 = 320 terms in the subsequent analysis.
The data used for comparison below "ANSYS" in all tables are computed by the finite element model in ANSYS Workbench (version 17.2).The middle surface model of the blade is established, and the finite element analysis is performed by using the SHELL181 linear quadrilateral element type.The layered section object in ANSYS mechanical is used to define the layered section of composite modeling.In order to get accurate results as much as possible, 20 layers in each skin (total of 41 layers in the thickness direction) are defined to simulate the cross section of the functionally graded sandwich blade.

Comparative Studies.
The homogeneous material blade is used first to demonstrate the numerical accuracy of the analytical model in this paper, and the first six nondimensional frequency parameters are compared with ANSYS and reference data [15] in Table 3.The results are in good agreement when the twist angle θ = 0 °to 45 °in the case of the aspect ratio l = 3.The results match well also when the twist angle θ = 0 °to 30 °in the case of the aspect ratio l = 1, but the maximum error (occurring at the 2nd order) is about 8.86% when θ increases to 45 °.The reason is that the Lame parameter A is approximated to 1 in order to improve the calculation speed in the process of numerical calculation.The error is smaller in the case of a larger aspect ratio and smaller twist angle, but the error becomes larger as the aspect ratio decreases and the twist angle increases.Thus, the aspect ratio l = 3 is used to obtain more accurate results.
The functionally graded sandwich blade is composed of metal and ceramic materials; their properties are given in Table 4.In Table 5, the first six nondimensional frequency parameters ω of the FG sandwich blade are compared with results in ANSYS Workbench, both showing good agreements, respectively.
The mode shapes corresponding to the first six orders of frequencies are also listed in Table 5.It should be noted that spanwise bending modes are divided into flapwise bending modes (whose displacements are predominantly in the direction normal to the middle surface) and edgewise bending modes (whose displacements are predominantly in the direction tangent to the surface) [15] in this paper.The modes 1B and 2B represent the first and second flapwise bending modes, respectively, and 1EB represents the first edgewise bending mode.The modes 1T and 2T represent the first and second torsional modes, respectively.The chordwise bending modes (CB, having two or more nodal lines approximately parallel to the pretwisted blade axis) are also found in the successive research.The curves such as 1B, 1T, and 2B and so on in Figures 6-8 are all obtained by using ANSYS Workbench.The curves such as 1st, 2nd, and 3rd are all obtained through numerical calculation.
The first six mode shapes of the FG sandwich blade with l = 3, h = 20, n = 2, h l = 1/2, and θ = 30 ∘ are given in Figure 9.It can be seen that the results obtained by the two methods are in good agreement.The results in Table 5 show that that the 1st modes are both of 1B mode shape when θ = 0 °, 15 °, and 30 °, but the 2nd mode is of 1T mode shape when θ = 0 °and 2B mode shape when θ = 15 °and 30 °.The same happens with the 1EB mode shape.The reason for this phenomenon is that the stiffness of the FG sandwich blade is changed due to the initial twist in (15).

4.3.
Free Vibrations of the FG Sandwich Blade.The effect of the twist angle on the vibration characteristics for the pretwist FG sandwich blade with geometric parameters as l = 3, h = 20, n = 2, and h l = 1/2 is studied.Table 6 gives the first six nondimensional frequency parameters ω of the FG sandwich blade with varying twist angles (θ from 0 °to 60 °), and the relationship between the nondimensional frequency parameters and twist angles for the FG sandwich blade is shown in Figure 6.It can be seen that the resulting curves are very close to the curves obtained by ANSYS in Figure 6(a).The first mode corresponds to the 1B vibration mode for all twist angles.The 2nd and 3rd nondimensional frequency parameter curves represent the 1T and 2B mode shapes when θ = 0 °, and they come close to each other and separate gradually with the twist angle increasing.The two curves are up to the nearest when the twist angle θ is about 5 °; the two modes exchange with each other, and the variation tendencies of nondimensional frequency parameters are also interchanged.
The variations for 4th to 6th can be found in Figure 6(b).
It should be noted that the FG sandwich blade designed by the geometric parameters corresponding to the frequency turning points will be prone to internal resonance.
The effect of the thickness ratio on the vibration characteristics for the pretwist FG sandwich blade is studied also.The first six nondimensional frequency parameters ω of the FG sandwich blade with varying thickness ratios (h from 10 to 100) is given in Table 7. Figure 7 shows the relationship between the nondimensional frequency parameters and the thickness ratio.It can be found that the nondimensional frequency parameters of 1B, 2B, 3B, and 4B modes remain essentially unchanged when the thickness ratio varies from 10 to 100.The nondimensional frequency parameters of 1T, 2T, and 1EB modes tend to increase with increasing thickness ratio.On the other side, the 3rd and 4th frequency parameter loci have obvious veering at h = 85, and the mode shapes exchanged between 1T and 3B.The same situation can be seen on the 5th and 6th modes, and higher-order frequency parameter locus veering happens more frequently.
The nondimensional frequency parameters ω for the first six modes of the FG sandwich blade with varying aspect ratios (l from 1 to 5) are given in Table 8 to study the effect of the aspect ratio on the vibration characteristics for the pretwist FG sandwich blade.Figure 8 shows the relationships between the first six frequency parameters, the corresponding mode shapes, and the aspect ratio.It is obvious that the 1st mode corresponds to the 1B vibration mode for all aspect ratios, and the value of the frequency parameter varies little.The 2nd and 3rd modes exchange with each other, and the   14 International Journal of Aerospace Engineering vibration tendencies of their frequency loci also exchange near l = 2.The change of the 4th to 6th nondimensional frequency parameters is relatively complicated.The mode shapes include 3B, 2T, 1EB, and 1CB.For example, the 4th to 6th mode shapes are 1CB, 2T, and 1EB when l = 1; 2T, 1EB, and 1CB when l = 1 5; and 1EB, 2T, and 3B when l = 1 75.All in all, the mode shape corresponding to the same order frequency is very sensitive to the variation of small aspect ratios, especially when l < 2. The tendencies of their frequency parameters become relatively stable until l > 2 5. Table 9 and Figure 10 show the variations of the first six nondimensional frequency parameters of the FG sandwich blade with the layer thickness ratio.Table 10 and Figure 11 show the relationships between the first six nondimensional frequency parameters and the scalar parameter of the volume fraction.It can be found that the content of ceramic increases with increasing layer thickness ratio and scalar parameter of the volume fraction, and the frequency parameters of all modes increase, respectively.
Table 11 shows the first six nondimensional frequency parameters of the FG sandwich blade for the nondimensional rotation velocity with the range of Ω 0 to 11 (about 15,000 rpm).A typical Campbell diagram is plotted to study the variation trend of the nondimensional frequency parameters for varying rotation velocities in Figure 12.Frequency locus veering is observed as expected: the 1st and 2nd frequency modes exchange between 1B and 1T around the rotation velocity Ω = 11; the 2nd and 3rd frequency modes exchange between 2B and 1T around Ω = 2; the 5th and 6th frequency modes exchange between 1EB and 2T around Ω = 3; and the 7th and 8th frequency modes exchange between 4B and 3T around Ω = 1 and exchange between 3T and 1CB again around the rotation velocity Ω = 9.In addition, the 8th frequency locus can be seen veering at Ω = 4 and the mode shape is changed from 4B to 1EB.It can be seen clearly that higher-order frequencies tend to have multiple frequency veering with the mode shapes exchanging and coupling.
Figure 13 shows the effects of the stagger angle for the first six nondimensional frequency parameters with varying rotation velocities.It can be found that the stagger angle has a predominant influence in the 1st and 2nd modes, and the effect is very marginal in higher modes.The effect of the stagger angle has little influence on all frequencies when the rotation velocity is low.The reason is that the centrifugal force has little effect on the stiffening and softening when the rotating speed is low, and the influence becomes more and more obvious with the rotating speed increasing.

Conclusion
In this paper, a dynamic model based on the thick shell theory is developed to investigate the free vibration behavior of functionally graded sandwich pretwisted blades.The validation of the homogeneous material and FG sandwich blade is performed by comparison to the literature and ANSYS results, both showing good agreement.
Interesting frequency veering and the associated mode shift phenomena are found and discussed for different twist angles, thickness ratios, and aspect ratios in static natural frequency analysis.The results show that the fundamental mode is flapwise bending in every case.The mode shapes of the 2nd to 6th orders are changing when different geometric parameters are used; especially, the thickness ratio is 15-50    International Journal of Aerospace Engineering and the aspect ratio is 1.0-2.5.Higher-order frequency loci (≥4) may veer more than once between the bending and torsional vibrations of the FG sandwich blade.The stagger angle is considered to study the dynamic frequency with varying rotation velocities.The results show that the frequency modes have more complicated shapes and stronger coupling in the rotating FG sandwich blade system.The Campbell diagram reveals that the rotating velocity, especially high rotating velocity, has a significant impact on the variations of frequency loci and mode shapes.

Figure 2 :
Figure 2: The dynamic coordinate system of the pretwisted rotating blade.

3
International Journal of Aerospace Engineering surface in this paper.The volume fraction of the metal phase is obtained as[27]

Figure 3 :
Figure 3: The section of the FG sandwich blade.

Figure 4 :
Figure 4: Effect of the scalar parameter of volume fraction in the FG sandwich blade.

Table 1 :
Convergence of nondimensional frequency parameters for the first six modes of the homogeneous material blade (l = 3, h = 20, and θ = 15 °).

Table 3 :
) with the initial twist angle.It can be found that both the nondimensional stiffness coefficients increase with increase in the initial twist angle.2.3.The Potential Energy of CentrifugalForce.The coordinate system transformation relationship between the xyz-coordinate and XYZ-coordinate is given by Comparison of nondimensional frequency parameters for the first six modes of homogeneous material blade (h = 20, μ = 0 3).

Table 4 :
The materials and properties of the FG sandwich blade.

Table 5 :
Comparisons of nondimensional frequency parameters for the first six modes of the FG sandwich blade (l = 3, h = 20, n = 2, and h l = 1/2).

Table 10 :
Nondimensional frequency parameters for the first six modes of the FG sandwich blade with varying scalar parameters of volume fraction (l = 3, h = 20, h l = 1/2, and θ = 30 °).