# Dynamic Analysis of Tapered Thin-Walled Beams Using Spectral Finite Element Method

Authors: Yiping Shen; Zhijun Zhu; Songlai Wang; Gang Wang
Journal: Shock and Vibration
 (2019)

---

## Abstract

Tapered thin-walled structures have been widely used in wind turbine and rotor blade. In this paper, a spectral finite element model is developed to investigate tapered thin-walled beam structures, in which torsion related warping effect is included. First, a set of fully coupled governing equations are derived using Hamilton’s principle to account for axial, bending, and torsion motion. Then, the differential transform method (DTM) is applied to obtain the semianalytical solutions in order to formulate the spectral finite element. Finally, numerical simulations are conducted for tapered thin-walled wind turbine rotor blades and validated by the ANSYS. Modal frequency results agree well with the ANSYS predictions, in which approximate 30,000 shell elements were used. In the SFEM, one single spectral finite element is needed to perform such calculations because the interpolation functions are deduced from the exact semianalytical solutions. Coupled axial-bending-torsion mode shapes are obtained as well. In summary, the proposed spectral finite element model is able to accurately and efficiently to perform the modal analysis for tapered thin-walled rotor blades. These modal frequency and mode shape results are important to carry out design and performance evaluation of the tapered thin-walled structures.


---

## Body

## 1. Introduction

Tapered thin-walled structures have been widely used in wind turbine, helicopter, and fix-wing aircraft due to their weight-saving and aerodynamic load carrying capabilities [1]. Accurate prediction of their dynamic behaviors is crucial to carry out structural design analysis and performance evaluation. Such structures are relatively long and slender, which are typically represented as thin-walled beams with complex closed cross section configurations. Normally, the shear center does not coincide with the centroid location in these thin-walled beams. Therefore, the torsion related warping effect cannot be neglected, which leads to a fully coupled axial-bending-torsion system. A tapered design is introduced to provide optimized load carrying capability and weight saving. However, due to the presence of these tapers, the governing differential equations contain variable coefficients. It is not possible to obtain closed-form solutions for the displacement field. Numerical approaches (e.g., finite element method) must be used to estimate the performance of tapered thin-walled rotor blade. Variable cross-sectional properties add challenge to characterize the behavior of these tapered thin-walled beam structures. The need of one-dimensional (1D) beam finite element model for tapered thin-walled beam still remains as an attractive research topic.The 1D beam model of the thin-walled structure is extensively established by Timoshenko [2], Vlasov [3], Gjelsvik [4], and Librescu and Song [1]. Many analyses have been formulated for symmetric and antisymmetric closed-section thin-walled beams with varying levels of assumptions, such as for axial warping and shell wall bending. Chandra [5] discussed the bending-twist coupling effects for box beams under bending, torsional, and extensional loads. Simo [6] extended Vlasov’s thin-walled beam theory to a fully nonlinear geometrically exact beam model. Lee [7] applied the same model to develop 1D model for flexural-torsional behavior of a laminated composite beam with an I-section. Ruta [8] proposed a direct 1D continuum beam model with respect to the centroidal axis and the shear center axis to describe the flexural-torsional buckling, and Pignataro and Rizzi [9] investigated the effects of warping on its buckling behavior. Librescu and Song’s work [1] contains the detailed formulations for 1D beam modeling of thin-walled structure, including axial, bending, torsion, and torsional-related warping, as well as the anisotropic thin-walled beam modeling of aircraft wing. Recently, improved 1D beam models are developed for thin-walled beam by Vo and Lee [10–12], Cárdenas [13], Zhang [14], and many others. In this paper, Vo and Lee’s model [10] is extended to model the tapered thin-walled blade.Based on aforementioned thin-walled beam models, corresponding conventional finite element methods (CFEMs) can be developed. For thin-walled rotor blades, 1D beam finite element is preferred compared to 2D shell element due to its efficiency. Many elements must be employed to discretize a thin-walled beam in order to account for the tapered configuration, because each element has different cross-sectional properties. Various efforts have been conducted to improve the displacement interpolation functions for tapered thin-walled beams [15–20]. Shin [21] presented a new tapered high-order element for tapered thin-walled box beams. Zappino [22] applied the Carrera unified formulation [23] to analyze tapered thin-walled composite structures.The spectral finite element method (SFEM) is a viable structural analysis approach that can provide high-fidelity predictions using comparatively small number of elements. It was initially introduced by Narayanan and Beskos [24] and also called the dynamic stiffness method [25, 26]. The SFEM is developed in the frequency domain [27, 28]. The discrete Fourier transform is applied to transform the governing equations from the time domain to the frequency domain. Natural frequencies are determined by solving transcendental equations instead of obtaining the eigenvalue problem solutions as practiced in the CFEM. Giurgiutiu and Stafford [29] solved the uncoupled equations for lag, flap, and pitch vibration of the uniform rotor blades by the Frobenius method. Huang [30] proposed a dynamic stiffness method for a rotating beam of nonuniform cross section by using the power series method. Wang [31] and Banerjee [32] applied the Frobenius method to obtain the semianalytical solutions for linearly tapered rotating beams under bending vibration and formulated the spectral finite element. Vinod [33] applied the exact solution of the transverse governing equation to develop an approximate spectral element for rotating Euler–Bernoulli beams in the frequency domain.The differential transform method (DTM) has been applied to solve the dynamic equation of a double tapered Timoshenko beam considering flapwise bending-torsion coupling recently [34]. The DTM is based on the Taylor series expansion, first introduced by Zhou in 1986 [35], which is an efficient method to solve such fully coupled governing equation with variable coefficients.In this paper, a spectral finite element model is developed to investigate tapered thin-walled rotor blade structures, in which torsion related warping effect is included. Vo and Lee’s model [11] is extended to formulate the kinematic equations. A set of fully coupled governing equations are derived using Hamilton’s principle to account for axial, bending, and torsion motion. Then, the differential transform method (DTM) is applied to obtain the semianalytical solutions in order to formulate the spectral finite element. Finally, modal frequency and mode shape results are obtained for tapered thin-walled wind turbine rotor blades and validated by the ANSYS.The paper is organized as follows. Section2 details the derivation of the axial-bending-torsion fully coupled governing equations, in which single-celled section is assumed and torsion-related warping effects are included. In Section 3, a SFEM is formulated for tapered rotor blade, in which the DTM is applied to obtain the semianalytical solutions. In Section 4, a tapered cylinder beam and a tapered wind turbine rotor blade are simulated to validate by ANSYS prediction. A summary and some concluding remarks are given in Section 5.


## 2. Modeling of Tapered Thin-Walled Rotor Blades

### 2.1. Kinematics

A generic single-cell tapered thin-walled rotor blade is shown in Figure1, in which the airfoil shape is DU40-A17 of a NREL 5 MW wind turbine [36]. The detailed assumptions and kinematic modeling for such thin-walled beam can be found in the works of Librescu and Song [1], Vo and Lee [11], and Zhang [15]; the development of such beam herein is presented in brief. Two sets of coordinate systems are used, as shown in Figures 1(a) and 1(b).Figure 1
(a) Geometry and coordinate system of a tapered thin-walled rotor blade. (b) Sectional coordinates and displacements.
(a)
(b)The first set is a global Cartesian coordinate systemoxyz. The z axis is the axial direction of the rotor blade, and its position is placed through the shear centers of the cross sections herein. The x and y axes are the transverse coordinates of the closed cross section. The second set is a local shell wall coordinate system zsn, where s is a contour coordinate measured along the middle surface of the shell wall and n is normal to the contour coordinate. Due to asymmetric airfoil configuration, the shear center does not coincide with the centroid of the thin-walled beam. Therefore, both flapwise and edgewise bending is coupled with torsion motion. Moreover, the torsion-related warping has a significant effect on the displacements and forces in axial, bending, and torsion motion. Key assumptions made in our model are listed as follows:(1)
The contour does not deform in its own plane(2)
The shear strain of the middle surface of the shell wall is to have the same distribution in the contour direction as in St.Venant torsion(3)
The Kirchhoff–Love assumption is applied to the thin wallAssumption 1 implies that the cross sections of the blade are assumed rigid in their own planes but are allowed to warp in the axial direction. Experiments and theoretical investigations have confirmed its essential accuracy for engineering applications [1, 4]. According to assumption 1, the middle surface displacements u¯ and v¯ at any point P in the local shell wall coordinate system can be expressed in terms of displacements U, V in the x and y directions, respectively, and the rotation angle ϕ about z axis in the global coordinate system [11, 15]:(1)u¯s,z=Uzsinθs−Vzcosθs−qsϕz,(2)v¯s,z=Uzcosθs+Vzsinθs+rsϕz.As shown in Figure1(b), the associated parameters rs and qs are given by(3)qs=xscosθs+yssinθs,rs=xssinθs−yscosθs,where θs is the angle between the x axis and the local s axis at the point s on the perimeter.The corresponding out-of-plane shell displacementw¯ can be obtained based on assumption 2. Neglecting the transverse shearing effects on the deformation of the blade, the shear strain γ¯sz of the middle surface can be expressed as(4)γ¯sz=∂v¯s,z∂z+∂w¯s,z∂s=ϕ′z∮rsds∮ds.Equation (4) gives the solution of the shell wall midsurface warping displacement. Substituting equation (2) into equation (4) and conducting integration of equation (4) with respect to s, the axial displacement can be determined as shown below:(5)w¯s,z=Wz−U′zxs−V′zys−ωsϕ′z,ωs=∫0srs−∮rsds/∮dsds,where dx=dssinθs and dy=dscosθs.According to the Kirchhoff–Love shell theory in assumption 3, the local displacements at any pointP on the shell wall can be expressed as(6)us,z,n=u¯s,z,vs,z,n=v¯s,z−n∂u¯s,z∂s,ws,z,n=w¯s,z−n∂u¯s,z∂z.The strains in the shell wall are given by(7)εz=ε¯z+nκ¯z,γsz=γ¯sz+nκ¯sz.(8)ε¯z=∂w¯s,z∂z=W′z−U″zxs−V″zys−ωsϕ″z,κ¯z=−∂2u¯s,z∂z2=−U″zsinθs+V″zcosθs+qsϕ″z,κ¯sz=−2∂2u¯s,z∂s∂z=2ϕ′z.Finally, the strains in the shell wall can be obtained by substituting equation (8) into equation (7) as shown below:(9)εz=W′z−xs+nsinθsU″z−ys−ncosθsV″z−ωs−nqsϕ″z,γsz=∮rsds∮ds+2nϕ′z.


### 2.2. Governing Equations

The total strain energyK is given by(10)K=12∫0L∮∫−t/2t/2Eεz2+Gγsz2dndsdz=12∫0L∮∫−t/2t/2σzεz+σszγszdndsdz.Substituting equation (9) into equation (10), the variation of strain energy can be rewritten in terms of general force and displacement components(11)δK=∫0LNzδW′−MyδU″−MxδV″−Mωδϕ″+Mtδϕ′dz,where Nz,My,Mx,Mt,andMω are the axial force, bending moments, torsional moment, and warping bimoment, respectively. They are defined by integrating over the cross section as(12)Nz=∮∫−t/2t/2σzdnds,My=∮∫−t/2t/2σzxs+nsinθdnds,Mx=∮∫−t/2t/2σzys−ncosθdnds,Mω=∮∫−t/2t/2σzωs−nqsdnds,Mt=∮∫−t/2t/2σsz∮rsds∮ds+2ndnds.The final expression of these force components are given below:(13)Nz=EAW′−SyU″−SxV″−Sωϕ″,My=ESyW′−IyU″−IxyV″−Iyωϕ″,Mx=ESxW′−IxyU″−IxV″−Ixωϕ″,Mω=ESωW′−IyωU″−IxωV″−Iωϕ″,Mt=GJϕ′.Corresponding sectional properties are varied alongz axis for a tapered thin-walled beam. These parameters at an arbitrary cross section can be simplified for a thin-walled beam and are shown below:(14)A=∮tds,Sy=∮xstds,Sx=∮ystds,Sω=∮ωstds,J=∮t∮rsds2∮ds2ds,Iy=∮x2stds,Ix=∮y2stds,Ixy=∮xsystds,Iω=∮ω2stds,Ixω=∮ysωstds,Iyω=∮xsωstds.It has been discussed by Gjelsvik [4] that the following parameters are zero when the pole O of the local sectional coordinate system is identical to the shear center:(15)Sω=Iyω=Ixω=0.There are two definitions of the shear center as reviewed in [37], named as kinematic and energetic shear centers. In essence, the kinematic shear center is used in this paper. The other parameters are calculated by using Shi’s finite segment method [38].The kinematic energy is given by(16)T=12∫0L∮∫−t/2t/2ρw˙2+u˙2+v˙2dndsdz=12∫0L∮∫−t/2t/2ρW˙−x+nsinθU˙′−y−ncosθV′˙−ω−nqϕ′˙2+U˙sinθ−V˙cosθ−qϕ˙2+U˙cosθ+V˙sinθ+rϕ˙2dndsdz.The variation of the kinematic energy is shown below:(17)δT=ρ∫0LAW˙−SyU˙′−SxV˙′δW˙+AU˙−Sxϕ˙δU˙+−SyW˙+IyU˙′+IxyV˙′δU˙′+AV˙+Syϕ˙δV˙+−SxW˙+IxyU˙′+IxV˙′δV˙′+Ix+Iyϕ˙−SxU˙+SyV˙δϕ˙+Iωϕ˙′δϕ˙′dz.Applying Hamilton’s principle, the governing equations and associated boundary conditions are obtained. The coupled axial-bending-torsion equations are(18)ρAW¨−ρSyU¨′−ρSxV¨′−EAW′′+ESyU″′+ESxV″′=0,(19)ρAU¨−ρSxϕ¨+ρSyW¨′−ρIyU¨′′−ρIxyV¨′′−ESyW′″+EIyU″″+EIxyV″″=0,(20)ρAV¨+ρSyϕ¨+ρSxW¨′−ρIxyU¨′′−ρIxV¨′′−ESxW′″+EIxyU″″+EIxV″″=0,(21)ρIx+Iyϕ¨−ρSxU¨+ρSyV¨−ρIωϕ¨′′−GJϕ′′+EIωϕ″″=0.The blade is a cantilever beam, corresponding boundary conditions are given by(22)z=0:W0=0,U0=0,U′0=0,V0=0,V′0=0,ϕ0=0,ϕ′0=0.(23)z=L:NzL=0,MyL=0,My′L=0,MxL=0,Mx′L=0,MtL−Mω′L=0,MωL=0.In this study, it is assumed that the chord lengthμz of rotor blade is uniformly tapered along z axis. It is defined by(24)μz=μ01+αz,where α is the taper ratio, α=μL/μ0−1/L and μL and μ0 are the chord length of its two ends, respectively.The sectional properties at any cross section can be rewritten as follows in order to account for the taper:(25)Az=A01+αz,Syz=Sy01+αz2,Sxz=Sx01+αz2,Jz=J01+αz3,Iyz=Iy01+αz3,Ixz=Ix01+αz3,Ixyz=Ixy01+αz3,Iωz=Iω01+αz5.Finally, the governing equations (18)–(21) can be modified to account for the variation of sectional properties:(26)ρA01+αzW¨−ρSy01+αz2U¨′−ρSx01+αz2V¨′−EA01+αzW′′+ESy01+αz2U″′+ESx01+αz2V″′=0,(27)ρA01+αzU¨−ρSx01+αz2ϕ¨+ρSy01+αz2W¨′−ρIy01+αz3U¨′′−ρIxy01+αz3V¨′′−ESy01+αz2W′″+EIy01+αz3U″″+EIxy01+αz3V″″=0,(28)ρA01+αzV¨+ρSy01+αz2ϕ¨+ρSx01+αz2W¨′−ρIxy01+αz3U¨′′−ρIx01+αz3V¨′′−ESx01+αz2W′″+EIxy01+αz3U″″+EIx01+αz3V″″=0,(29)ρIx0+Iy01+αz3ϕ¨−ρSx01+αz2U¨+ρSy01+αz2V¨−ρIω01+αz5ϕ¨′′−GJ01+αz3ϕ′′+EIω01+αz5ϕ″″=0.


## 2.1. Kinematics

A generic single-cell tapered thin-walled rotor blade is shown in Figure1, in which the airfoil shape is DU40-A17 of a NREL 5 MW wind turbine [36]. The detailed assumptions and kinematic modeling for such thin-walled beam can be found in the works of Librescu and Song [1], Vo and Lee [11], and Zhang [15]; the development of such beam herein is presented in brief. Two sets of coordinate systems are used, as shown in Figures 1(a) and 1(b).Figure 1
(a) Geometry and coordinate system of a tapered thin-walled rotor blade. (b) Sectional coordinates and displacements.
(a)
(b)The first set is a global Cartesian coordinate systemoxyz. The z axis is the axial direction of the rotor blade, and its position is placed through the shear centers of the cross sections herein. The x and y axes are the transverse coordinates of the closed cross section. The second set is a local shell wall coordinate system zsn, where s is a contour coordinate measured along the middle surface of the shell wall and n is normal to the contour coordinate. Due to asymmetric airfoil configuration, the shear center does not coincide with the centroid of the thin-walled beam. Therefore, both flapwise and edgewise bending is coupled with torsion motion. Moreover, the torsion-related warping has a significant effect on the displacements and forces in axial, bending, and torsion motion. Key assumptions made in our model are listed as follows:(1)
The contour does not deform in its own plane(2)
The shear strain of the middle surface of the shell wall is to have the same distribution in the contour direction as in St.Venant torsion(3)
The Kirchhoff–Love assumption is applied to the thin wallAssumption 1 implies that the cross sections of the blade are assumed rigid in their own planes but are allowed to warp in the axial direction. Experiments and theoretical investigations have confirmed its essential accuracy for engineering applications [1, 4]. According to assumption 1, the middle surface displacements u¯ and v¯ at any point P in the local shell wall coordinate system can be expressed in terms of displacements U, V in the x and y directions, respectively, and the rotation angle ϕ about z axis in the global coordinate system [11, 15]:(1)u¯s,z=Uzsinθs−Vzcosθs−qsϕz,(2)v¯s,z=Uzcosθs+Vzsinθs+rsϕz.As shown in Figure1(b), the associated parameters rs and qs are given by(3)qs=xscosθs+yssinθs,rs=xssinθs−yscosθs,where θs is the angle between the x axis and the local s axis at the point s on the perimeter.The corresponding out-of-plane shell displacementw¯ can be obtained based on assumption 2. Neglecting the transverse shearing effects on the deformation of the blade, the shear strain γ¯sz of the middle surface can be expressed as(4)γ¯sz=∂v¯s,z∂z+∂w¯s,z∂s=ϕ′z∮rsds∮ds.Equation (4) gives the solution of the shell wall midsurface warping displacement. Substituting equation (2) into equation (4) and conducting integration of equation (4) with respect to s, the axial displacement can be determined as shown below:(5)w¯s,z=Wz−U′zxs−V′zys−ωsϕ′z,ωs=∫0srs−∮rsds/∮dsds,where dx=dssinθs and dy=dscosθs.According to the Kirchhoff–Love shell theory in assumption 3, the local displacements at any pointP on the shell wall can be expressed as(6)us,z,n=u¯s,z,vs,z,n=v¯s,z−n∂u¯s,z∂s,ws,z,n=w¯s,z−n∂u¯s,z∂z.The strains in the shell wall are given by(7)εz=ε¯z+nκ¯z,γsz=γ¯sz+nκ¯sz.(8)ε¯z=∂w¯s,z∂z=W′z−U″zxs−V″zys−ωsϕ″z,κ¯z=−∂2u¯s,z∂z2=−U″zsinθs+V″zcosθs+qsϕ″z,κ¯sz=−2∂2u¯s,z∂s∂z=2ϕ′z.Finally, the strains in the shell wall can be obtained by substituting equation (8) into equation (7) as shown below:(9)εz=W′z−xs+nsinθsU″z−ys−ncosθsV″z−ωs−nqsϕ″z,γsz=∮rsds∮ds+2nϕ′z.


## 2.2. Governing Equations

The total strain energyK is given by(10)K=12∫0L∮∫−t/2t/2Eεz2+Gγsz2dndsdz=12∫0L∮∫−t/2t/2σzεz+σszγszdndsdz.Substituting equation (9) into equation (10), the variation of strain energy can be rewritten in terms of general force and displacement components(11)δK=∫0LNzδW′−MyδU″−MxδV″−Mωδϕ″+Mtδϕ′dz,where Nz,My,Mx,Mt,andMω are the axial force, bending moments, torsional moment, and warping bimoment, respectively. They are defined by integrating over the cross section as(12)Nz=∮∫−t/2t/2σzdnds,My=∮∫−t/2t/2σzxs+nsinθdnds,Mx=∮∫−t/2t/2σzys−ncosθdnds,Mω=∮∫−t/2t/2σzωs−nqsdnds,Mt=∮∫−t/2t/2σsz∮rsds∮ds+2ndnds.The final expression of these force components are given below:(13)Nz=EAW′−SyU″−SxV″−Sωϕ″,My=ESyW′−IyU″−IxyV″−Iyωϕ″,Mx=ESxW′−IxyU″−IxV″−Ixωϕ″,Mω=ESωW′−IyωU″−IxωV″−Iωϕ″,Mt=GJϕ′.Corresponding sectional properties are varied alongz axis for a tapered thin-walled beam. These parameters at an arbitrary cross section can be simplified for a thin-walled beam and are shown below:(14)A=∮tds,Sy=∮xstds,Sx=∮ystds,Sω=∮ωstds,J=∮t∮rsds2∮ds2ds,Iy=∮x2stds,Ix=∮y2stds,Ixy=∮xsystds,Iω=∮ω2stds,Ixω=∮ysωstds,Iyω=∮xsωstds.It has been discussed by Gjelsvik [4] that the following parameters are zero when the pole O of the local sectional coordinate system is identical to the shear center:(15)Sω=Iyω=Ixω=0.There are two definitions of the shear center as reviewed in [37], named as kinematic and energetic shear centers. In essence, the kinematic shear center is used in this paper. The other parameters are calculated by using Shi’s finite segment method [38].The kinematic energy is given by(16)T=12∫0L∮∫−t/2t/2ρw˙2+u˙2+v˙2dndsdz=12∫0L∮∫−t/2t/2ρW˙−x+nsinθU˙′−y−ncosθV′˙−ω−nqϕ′˙2+U˙sinθ−V˙cosθ−qϕ˙2+U˙cosθ+V˙sinθ+rϕ˙2dndsdz.The variation of the kinematic energy is shown below:(17)δT=ρ∫0LAW˙−SyU˙′−SxV˙′δW˙+AU˙−Sxϕ˙δU˙+−SyW˙+IyU˙′+IxyV˙′δU˙′+AV˙+Syϕ˙δV˙+−SxW˙+IxyU˙′+IxV˙′δV˙′+Ix+Iyϕ˙−SxU˙+SyV˙δϕ˙+Iωϕ˙′δϕ˙′dz.Applying Hamilton’s principle, the governing equations and associated boundary conditions are obtained. The coupled axial-bending-torsion equations are(18)ρAW¨−ρSyU¨′−ρSxV¨′−EAW′′+ESyU″′+ESxV″′=0,(19)ρAU¨−ρSxϕ¨+ρSyW¨′−ρIyU¨′′−ρIxyV¨′′−ESyW′″+EIyU″″+EIxyV″″=0,(20)ρAV¨+ρSyϕ¨+ρSxW¨′−ρIxyU¨′′−ρIxV¨′′−ESxW′″+EIxyU″″+EIxV″″=0,(21)ρIx+Iyϕ¨−ρSxU¨+ρSyV¨−ρIωϕ¨′′−GJϕ′′+EIωϕ″″=0.The blade is a cantilever beam, corresponding boundary conditions are given by(22)z=0:W0=0,U0=0,U′0=0,V0=0,V′0=0,ϕ0=0,ϕ′0=0.(23)z=L:NzL=0,MyL=0,My′L=0,MxL=0,Mx′L=0,MtL−Mω′L=0,MωL=0.In this study, it is assumed that the chord lengthμz of rotor blade is uniformly tapered along z axis. It is defined by(24)μz=μ01+αz,where α is the taper ratio, α=μL/μ0−1/L and μL and μ0 are the chord length of its two ends, respectively.The sectional properties at any cross section can be rewritten as follows in order to account for the taper:(25)Az=A01+αz,Syz=Sy01+αz2,Sxz=Sx01+αz2,Jz=J01+αz3,Iyz=Iy01+αz3,Ixz=Ix01+αz3,Ixyz=Ixy01+αz3,Iωz=Iω01+αz5.Finally, the governing equations (18)–(21) can be modified to account for the variation of sectional properties:(26)ρA01+αzW¨−ρSy01+αz2U¨′−ρSx01+αz2V¨′−EA01+αzW′′+ESy01+αz2U″′+ESx01+αz2V″′=0,(27)ρA01+αzU¨−ρSx01+αz2ϕ¨+ρSy01+αz2W¨′−ρIy01+αz3U¨′′−ρIxy01+αz3V¨′′−ESy01+αz2W′″+EIy01+αz3U″″+EIxy01+αz3V″″=0,(28)ρA01+αzV¨+ρSy01+αz2ϕ¨+ρSx01+αz2W¨′−ρIxy01+αz3U¨′′−ρIx01+αz3V¨′′−ESx01+αz2W′″+EIxy01+αz3U″″+EIx01+αz3V″″=0,(29)ρIx0+Iy01+αz3ϕ¨−ρSx01+αz2U¨+ρSy01+αz2V¨−ρIω01+αz5ϕ¨′′−GJ01+αz3ϕ′′+EIω01+αz5ϕ″″=0.


## 3. Spectral Finite Element Model

In this section, we will introduce key steps to formulate the spectral finite element. First, semianalytical solutions of the governing equations are determined in the frequency domain using the DFM. Then, similar to the conventional finite element method (CFEM), the nodal displacement and force vector are introduced to a dynamic stiffness matrix which is derived to relate the nodal displacement vector with nodal force vector. Details will be presented in the following subsections.


### 3.1. Semianalytical Solutions

Assuming harmonic responses, the displacementsU,V,W,andϕ in the time domain can be written in the following forms [39]:(30)Uz,t=Uazejωnt,Vz,t=Vazejωnt,Wz,t=Wazejωnt,ϕz,t=ϕazejωnt,where Ua,Va,Wa,andϕa are the steady state frequency dependent displacement solutions. ωn is the angular excitation frequency.The governing equations (26)–(29) can be rewritten as(31)−ρA0ωn21+αzWa−EA0αWa′+1+αzWa″+ρSy0ωn21+αz2Ua′+ESy02α1+αzUa″+1+αz2Ua″′+ρSx0ωn21+αz2Va′+ESx02α1+αzVa″+1+αz2Va″′=0,(32)−ρSy0ωn22α1+αzWa+1+αz2Wa′−ESy02α2Wa′+4α1+αzWa″+1+αz2Wa″′−ρA0ωn21+αzUa+ρIy0ωn23α1+αz2Ua′+1+αz3Ua″+EIy06α21+αzUa″+6α1+αz2Ua″′+1+αz3Ua″′+ρIxy0ωn23α1+αz2Va′+1+αz3Va″+EIxy06α21+αzVa″+6α1+αz2Va″′+1+αz3Va″′+ρSx0ωn21+αz2ϕa=0,(33)−ρSx0ωn22α1+αzWa+1+αz2Wa′−ESx,02α2Wa′+4α1+αzWa″+1+αz2Wa″′+ρIxy0ωn23α1+αz2Ua′+1+αz3Ua″+EIxy06α21+αzUa″+6α1+αz2Ua″′+1+αz3Ua″′−ρA0ωn21+αzVa+ρIx0ωn23α1+αz2Va′+1+αz3Va″+EIx06α21+αzVa″+6α1+αz2Va″′+1+αz3Va″′−ρSy0ωn21+αz2ϕa=0,(34)ρSx0ωn2Ua−ρSy0ωn2Va−ρωn2Ix0+Iy01+αzϕa+ρIω0ωn25α1+αz2ϕa′+1+αz3ϕa″−GJ03αϕa′+1+αzϕa″+EIω020α21+αzϕa″+10α1+αz2ϕa″′+1+αz3ϕa″″=0.Note thatWa″′ can be obtained by calculating the differential of Wa″ from equation (31). It yields(35)Wa″′=−ρωn2EWa′+α21+αz2Wa′−α1+αzWa″+ρSy0ωn2EA0αUa′+1+αzUa″+Sy0A03αUa″′+1+αzUa″″+ρSx0ωn2EA0αVa′+1+αzVa″+Sx0A03αVa″′+1+αzVa″″.Substituting equation (35) into equations (32) and (33), the corresponding governing equations can be rewritten as(36)−ρSy0ωn22α1+αzWa−ESy03α2Wa′+3α1+αzWa″−ρA0ωn21+αzUa+ρωn2α3Iy0−S˜y01+αz2Ua′+Iy0−S˜y01+αz3Ua″+E6α2Iy01+αzUa″+3α2Iy0−S˜y01+αz2Ua″′+Iy0−S˜y01+αz3Ua″″+ρωn2α3Ixy0−S˜xy01+αz2Va′+Ixy0−S˜xy01+αz3Va″+E6α2Ixy01+αzVa″+3α2Ixy0−S˜xy01+αz2Va″′+Ixy0−S˜xy01+αz3Va″″+ρSx0ωn21+αz2ϕa=0,(37)−ρSx0ωn22α1+αzWa−ESx,03α2Wa′+3α1+αzWa″+ρωn2α3Ixy0−S˜xy01+αz2Ua′+Ixy0−S˜xy01+αz3Ua″+E6α2Ixy01+αzUa″+3α2Ixy0−S˜xy01+αz2Ua″′+Ixy0−S˜xy01+αz3Ua″″−ρA0ωn21+αzVa+ρωn2α3Ix0−S˜x01+αz2Va′+Ix0−S˜x01+αz3Va″+E6α2Ix01+αzVa″+3α2Ix0−S˜x01+αz2Va″′+Ix0−S˜x01+αz3Va″″−ρSy0ωn21+αz2ϕa=0,where S˜x0=Sx02/A0,S˜y0=Sy02/A0,S˜xy0=Sx0Sy0/A0.


#### 3.1.1. Definition in the DFM

The DFM is an efficient method to obtain semianalytical solution of the governing differential equations. In this method, certain transformation rules are applied; the governing equations and the boundary conditions are transformed into a set of algebraic equations, which is an iterative procedure to obtain analytic Taylor Series solutions of the governing differential equations. The basic definitions and the application procedure of DFM are presented in the following.In DFM, the differential transform of the functionfz is defined as(38)Fk=1k!dkfzdzkz=z0,where fz is the original function and Fk is the transformed function.The inverse transformation can be expressed as(39)fz=∑k=0pFkz−z0k=∑k=0pz−z0kk!dkfzdzkz=z0,where p denotes the finite terms Fk of the transformed function, which depends on the convergence rate of the modal frequency prediction of rotor blade.Theorems that are used in the transformation of the governing differential equations and the boundary conditions of rotor blade are expressed as [34](40)fz=gz±hz⟹Fk=Gk±Hk,fz=λgz⟹Fk=λGk,fz=gzhz⟹Fk=∑m=0kGk−mHm,fz=dngzdzn⟹Fk=k+n!k!Gk+n,fz=zn⟹Fk=δk,m=0,ifk≠m,1,ifk=m.


#### 3.1.2. Transformations and Solutions of the Governing Equations

Based on the theorems in equation (40), the governing equations (31), (34), (36), and (37) can be rewritten in an iterative format. Because of the limited space, only the transformation of equation (31) is given in the following:(41)−ρA0ωn2∑m=0kδm,0+αδm,1Wak−m−EA0k+1!k!αWak+1+k+2!k!∑m=0kδm,0+αδm,1Wak−m+2+ρSy0ωn2k+1!k!∑m=0kδm,0+2αδm,1+α2δm,2Uak−m+1+ESy,02αk+2!k!∑m=0kδm,0+αδm,1Uak−m+2+k+3!k!∑m=0kδm,0+2αδm,1+α2δm,2Uak−m+3+ρSx0ωn2k+1!k!∑m=0kδm,0+2αδm,1+α2δm,2Vak−m+1+ESx,02αk+2!k!∑m=0kδm,0+αδm,1Vak−m+2+k+3!k!∑m=0kδm,0+2αδm,1+α2δm,2Vak−m+3=0.Furthermore, the boundary conditions in equation (22) can be rewritten as(42)Wa0=0,Ua0=0,Va0=0,ϕa0=0.Considering equation (13), the boundary conditions shown in equation (24) are yielded as the following equations:(43)A0∑k=0pk+1!k!Wak+1Lk−Sy01+αL∑k=0pk+2!k!Uak+2Lk−Sx01+αL∑k=0pk+2!k!Vak+2Lk=0,(44)Sy0∑k=0pk+1!k!Wak+1Lk−Iy01+αL∑k=0pk+2!k!Uak+2Lk−Ixy01+αL∑k=0pk+2!k!Vak+2Lk=0,(45)Sy02α∑k=0pk+1!k!Wak+1Lk+1+αL∑k=0pk+2!k!Wak+2Lk−Iy03α1+αL∑k=0pk+2!k!Uak+2Lk+1+αL2∑k=0pk+3!k!Uak+3Lk−Ixy03α1+αL∑k=0pk+2!k!Vak+2Lk+1+αL2∑k=0pk+3!k!Vak+3Lk=0,(46)Sx0∑k=0pk+1!k!Wak+1Lk−Ixy01+αL∑k=0pk+2!k!Uak+2Lk−Ix01+αL∑k=0pk+2!k!Vak+2Lk=0,(47)Sx02α∑k=0pk+1!k!Wak+1Lk+1+αL∑k=0pk+2!k!Wak+2Lk−Ixy03α1+αL∑k=0pk+2!k!Uak+2Lk+1+αL2∑k=0pk+3!k!Uak+3Lk−Ix03α1+αL∑k=0pk+2!k!Vak+2Lk+1+αL2∑k=0pk+3!k!Vak+3Lk=0,(48)GJ0∑k=0pk+1!k!ϕak+1Lk−EIω05α1+αL∑k=0pk+2!k!ϕak+2Lk+1+αL2∑k=0pk+3!k!ϕak+3Lk=0,(49)∑k=0pk+2!k!ϕak+2Lk=0.AssumeWa1=c1,Ua2=c2,Ua3=c3,Va2=c4,Va3=c5,ϕa2=c6,ϕa3=c7 and substitute the boundary conditions expressed into the governing equations rewritten in the iterative format using DTM; one of them is presented in equation (41); the transformed function Wak, Uak, Vak, and ϕak can be obtained with iterative process. Consequently, exact solutions of Wa, Ua, Va, and ϕa can be expressed as the sum of all the terms of the corresponding transformed function. The natural frequencies and mode shapes can be determined by applying the boundary conditions and solving a transcendental equation to obtain the unknown ci, i=1,…,7, which is yielded as(50)H11pH12p⋯H17pH21pH22p⋯H27p⋯⋯⋯⋯H71pH72p⋯H77pc1c2⋯c7=0,where Hijpi,j=1,2…7 are polynomials about unknown frequency ωn.Therefore, the eigenvalues can be calculated by taking the determinant of the matrix to be equal to zero:(51)H11pH12p⋯H17pH21pH22p⋯H27p⋯⋯⋯⋯H71pH72p⋯H77p=0.Solving equation (51), the modal frequencies ωn can be obtained. It is necessary to take more terms of Wak, Uak, Vak, and ϕak to evaluate the low order mode frequencies with satisfied precision. This procedure using DFM is different from the CFEM; the corresponding shape function is semianalytical solutions of the governing equations. The iterative calculation for the terms of transformed function Wak, Uak, Vak, and ϕak and modal frequencies of rotor blade can be quickly solved with a symbolic computational software.


### 3.2. Spectral Finite Element Formulation

Similar to the conventional finite element, a dynamic stiffness matrix can be defined to relate the nodal displacements with nodal forces. Figure2 shows the tapered spectral finite element for a rotor blade, in which both nodal displacement and force vectors are defined in the frequency domain, denoted by u^ and F^, respectively(52)u^=U1V1W1ϕ1U2V2W2ϕ2,F^=My1Mx1N1Mt1My2Mx2N2Mt2.Figure 2
Spectral finite element of a tapered thin-walled rotor blade.The nodal force vector can be related with the nodal displacement using a dynamic stiffness matrix,D, which is defined as the following [39]:(53)F^=Dωnu^.Consequently, such standard matrix representation of elemental dynamic stiffness can be used to assemble elements as used in the CFEM. By applying appropriate boundary conditions, we can conduct dynamic analysis of taped thin-walled rotor blade.


## 3.1. Semianalytical Solutions

Assuming harmonic responses, the displacementsU,V,W,andϕ in the time domain can be written in the following forms [39]:(30)Uz,t=Uazejωnt,Vz,t=Vazejωnt,Wz,t=Wazejωnt,ϕz,t=ϕazejωnt,where Ua,Va,Wa,andϕa are the steady state frequency dependent displacement solutions. ωn is the angular excitation frequency.The governing equations (26)–(29) can be rewritten as(31)−ρA0ωn21+αzWa−EA0αWa′+1+αzWa″+ρSy0ωn21+αz2Ua′+ESy02α1+αzUa″+1+αz2Ua″′+ρSx0ωn21+αz2Va′+ESx02α1+αzVa″+1+αz2Va″′=0,(32)−ρSy0ωn22α1+αzWa+1+αz2Wa′−ESy02α2Wa′+4α1+αzWa″+1+αz2Wa″′−ρA0ωn21+αzUa+ρIy0ωn23α1+αz2Ua′+1+αz3Ua″+EIy06α21+αzUa″+6α1+αz2Ua″′+1+αz3Ua″′+ρIxy0ωn23α1+αz2Va′+1+αz3Va″+EIxy06α21+αzVa″+6α1+αz2Va″′+1+αz3Va″′+ρSx0ωn21+αz2ϕa=0,(33)−ρSx0ωn22α1+αzWa+1+αz2Wa′−ESx,02α2Wa′+4α1+αzWa″+1+αz2Wa″′+ρIxy0ωn23α1+αz2Ua′+1+αz3Ua″+EIxy06α21+αzUa″+6α1+αz2Ua″′+1+αz3Ua″′−ρA0ωn21+αzVa+ρIx0ωn23α1+αz2Va′+1+αz3Va″+EIx06α21+αzVa″+6α1+αz2Va″′+1+αz3Va″′−ρSy0ωn21+αz2ϕa=0,(34)ρSx0ωn2Ua−ρSy0ωn2Va−ρωn2Ix0+Iy01+αzϕa+ρIω0ωn25α1+αz2ϕa′+1+αz3ϕa″−GJ03αϕa′+1+αzϕa″+EIω020α21+αzϕa″+10α1+αz2ϕa″′+1+αz3ϕa″″=0.Note thatWa″′ can be obtained by calculating the differential of Wa″ from equation (31). It yields(35)Wa″′=−ρωn2EWa′+α21+αz2Wa′−α1+αzWa″+ρSy0ωn2EA0αUa′+1+αzUa″+Sy0A03αUa″′+1+αzUa″″+ρSx0ωn2EA0αVa′+1+αzVa″+Sx0A03αVa″′+1+αzVa″″.Substituting equation (35) into equations (32) and (33), the corresponding governing equations can be rewritten as(36)−ρSy0ωn22α1+αzWa−ESy03α2Wa′+3α1+αzWa″−ρA0ωn21+αzUa+ρωn2α3Iy0−S˜y01+αz2Ua′+Iy0−S˜y01+αz3Ua″+E6α2Iy01+αzUa″+3α2Iy0−S˜y01+αz2Ua″′+Iy0−S˜y01+αz3Ua″″+ρωn2α3Ixy0−S˜xy01+αz2Va′+Ixy0−S˜xy01+αz3Va″+E6α2Ixy01+αzVa″+3α2Ixy0−S˜xy01+αz2Va″′+Ixy0−S˜xy01+αz3Va″″+ρSx0ωn21+αz2ϕa=0,(37)−ρSx0ωn22α1+αzWa−ESx,03α2Wa′+3α1+αzWa″+ρωn2α3Ixy0−S˜xy01+αz2Ua′+Ixy0−S˜xy01+αz3Ua″+E6α2Ixy01+αzUa″+3α2Ixy0−S˜xy01+αz2Ua″′+Ixy0−S˜xy01+αz3Ua″″−ρA0ωn21+αzVa+ρωn2α3Ix0−S˜x01+αz2Va′+Ix0−S˜x01+αz3Va″+E6α2Ix01+αzVa″+3α2Ix0−S˜x01+αz2Va″′+Ix0−S˜x01+αz3Va″″−ρSy0ωn21+αz2ϕa=0,where S˜x0=Sx02/A0,S˜y0=Sy02/A0,S˜xy0=Sx0Sy0/A0.


### 3.1.1. Definition in the DFM

The DFM is an efficient method to obtain semianalytical solution of the governing differential equations. In this method, certain transformation rules are applied; the governing equations and the boundary conditions are transformed into a set of algebraic equations, which is an iterative procedure to obtain analytic Taylor Series solutions of the governing differential equations. The basic definitions and the application procedure of DFM are presented in the following.In DFM, the differential transform of the functionfz is defined as(38)Fk=1k!dkfzdzkz=z0,where fz is the original function and Fk is the transformed function.The inverse transformation can be expressed as(39)fz=∑k=0pFkz−z0k=∑k=0pz−z0kk!dkfzdzkz=z0,where p denotes the finite terms Fk of the transformed function, which depends on the convergence rate of the modal frequency prediction of rotor blade.Theorems that are used in the transformation of the governing differential equations and the boundary conditions of rotor blade are expressed as [34](40)fz=gz±hz⟹Fk=Gk±Hk,fz=λgz⟹Fk=λGk,fz=gzhz⟹Fk=∑m=0kGk−mHm,fz=dngzdzn⟹Fk=k+n!k!Gk+n,fz=zn⟹Fk=δk,m=0,ifk≠m,1,ifk=m.


### 3.1.2. Transformations and Solutions of the Governing Equations

Based on the theorems in equation (40), the governing equations (31), (34), (36), and (37) can be rewritten in an iterative format. Because of the limited space, only the transformation of equation (31) is given in the following:(41)−ρA0ωn2∑m=0kδm,0+αδm,1Wak−m−EA0k+1!k!αWak+1+k+2!k!∑m=0kδm,0+αδm,1Wak−m+2+ρSy0ωn2k+1!k!∑m=0kδm,0+2αδm,1+α2δm,2Uak−m+1+ESy,02αk+2!k!∑m=0kδm,0+αδm,1Uak−m+2+k+3!k!∑m=0kδm,0+2αδm,1+α2δm,2Uak−m+3+ρSx0ωn2k+1!k!∑m=0kδm,0+2αδm,1+α2δm,2Vak−m+1+ESx,02αk+2!k!∑m=0kδm,0+αδm,1Vak−m+2+k+3!k!∑m=0kδm,0+2αδm,1+α2δm,2Vak−m+3=0.Furthermore, the boundary conditions in equation (22) can be rewritten as(42)Wa0=0,Ua0=0,Va0=0,ϕa0=0.Considering equation (13), the boundary conditions shown in equation (24) are yielded as the following equations:(43)A0∑k=0pk+1!k!Wak+1Lk−Sy01+αL∑k=0pk+2!k!Uak+2Lk−Sx01+αL∑k=0pk+2!k!Vak+2Lk=0,(44)Sy0∑k=0pk+1!k!Wak+1Lk−Iy01+αL∑k=0pk+2!k!Uak+2Lk−Ixy01+αL∑k=0pk+2!k!Vak+2Lk=0,(45)Sy02α∑k=0pk+1!k!Wak+1Lk+1+αL∑k=0pk+2!k!Wak+2Lk−Iy03α1+αL∑k=0pk+2!k!Uak+2Lk+1+αL2∑k=0pk+3!k!Uak+3Lk−Ixy03α1+αL∑k=0pk+2!k!Vak+2Lk+1+αL2∑k=0pk+3!k!Vak+3Lk=0,(46)Sx0∑k=0pk+1!k!Wak+1Lk−Ixy01+αL∑k=0pk+2!k!Uak+2Lk−Ix01+αL∑k=0pk+2!k!Vak+2Lk=0,(47)Sx02α∑k=0pk+1!k!Wak+1Lk+1+αL∑k=0pk+2!k!Wak+2Lk−Ixy03α1+αL∑k=0pk+2!k!Uak+2Lk+1+αL2∑k=0pk+3!k!Uak+3Lk−Ix03α1+αL∑k=0pk+2!k!Vak+2Lk+1+αL2∑k=0pk+3!k!Vak+3Lk=0,(48)GJ0∑k=0pk+1!k!ϕak+1Lk−EIω05α1+αL∑k=0pk+2!k!ϕak+2Lk+1+αL2∑k=0pk+3!k!ϕak+3Lk=0,(49)∑k=0pk+2!k!ϕak+2Lk=0.AssumeWa1=c1,Ua2=c2,Ua3=c3,Va2=c4,Va3=c5,ϕa2=c6,ϕa3=c7 and substitute the boundary conditions expressed into the governing equations rewritten in the iterative format using DTM; one of them is presented in equation (41); the transformed function Wak, Uak, Vak, and ϕak can be obtained with iterative process. Consequently, exact solutions of Wa, Ua, Va, and ϕa can be expressed as the sum of all the terms of the corresponding transformed function. The natural frequencies and mode shapes can be determined by applying the boundary conditions and solving a transcendental equation to obtain the unknown ci, i=1,…,7, which is yielded as(50)H11pH12p⋯H17pH21pH22p⋯H27p⋯⋯⋯⋯H71pH72p⋯H77pc1c2⋯c7=0,where Hijpi,j=1,2…7 are polynomials about unknown frequency ωn.Therefore, the eigenvalues can be calculated by taking the determinant of the matrix to be equal to zero:(51)H11pH12p⋯H17pH21pH22p⋯H27p⋯⋯⋯⋯H71pH72p⋯H77p=0.Solving equation (51), the modal frequencies ωn can be obtained. It is necessary to take more terms of Wak, Uak, Vak, and ϕak to evaluate the low order mode frequencies with satisfied precision. This procedure using DFM is different from the CFEM; the corresponding shape function is semianalytical solutions of the governing equations. The iterative calculation for the terms of transformed function Wak, Uak, Vak, and ϕak and modal frequencies of rotor blade can be quickly solved with a symbolic computational software.


## 3.1.1. Definition in the DFM

The DFM is an efficient method to obtain semianalytical solution of the governing differential equations. In this method, certain transformation rules are applied; the governing equations and the boundary conditions are transformed into a set of algebraic equations, which is an iterative procedure to obtain analytic Taylor Series solutions of the governing differential equations. The basic definitions and the application procedure of DFM are presented in the following.In DFM, the differential transform of the functionfz is defined as(38)Fk=1k!dkfzdzkz=z0,where fz is the original function and Fk is the transformed function.The inverse transformation can be expressed as(39)fz=∑k=0pFkz−z0k=∑k=0pz−z0kk!dkfzdzkz=z0,where p denotes the finite terms Fk of the transformed function, which depends on the convergence rate of the modal frequency prediction of rotor blade.Theorems that are used in the transformation of the governing differential equations and the boundary conditions of rotor blade are expressed as [34](40)fz=gz±hz⟹Fk=Gk±Hk,fz=λgz⟹Fk=λGk,fz=gzhz⟹Fk=∑m=0kGk−mHm,fz=dngzdzn⟹Fk=k+n!k!Gk+n,fz=zn⟹Fk=δk,m=0,ifk≠m,1,ifk=m.


## 3.1.2. Transformations and Solutions of the Governing Equations

Based on the theorems in equation (40), the governing equations (31), (34), (36), and (37) can be rewritten in an iterative format. Because of the limited space, only the transformation of equation (31) is given in the following:(41)−ρA0ωn2∑m=0kδm,0+αδm,1Wak−m−EA0k+1!k!αWak+1+k+2!k!∑m=0kδm,0+αδm,1Wak−m+2+ρSy0ωn2k+1!k!∑m=0kδm,0+2αδm,1+α2δm,2Uak−m+1+ESy,02αk+2!k!∑m=0kδm,0+αδm,1Uak−m+2+k+3!k!∑m=0kδm,0+2αδm,1+α2δm,2Uak−m+3+ρSx0ωn2k+1!k!∑m=0kδm,0+2αδm,1+α2δm,2Vak−m+1+ESx,02αk+2!k!∑m=0kδm,0+αδm,1Vak−m+2+k+3!k!∑m=0kδm,0+2αδm,1+α2δm,2Vak−m+3=0.Furthermore, the boundary conditions in equation (22) can be rewritten as(42)Wa0=0,Ua0=0,Va0=0,ϕa0=0.Considering equation (13), the boundary conditions shown in equation (24) are yielded as the following equations:(43)A0∑k=0pk+1!k!Wak+1Lk−Sy01+αL∑k=0pk+2!k!Uak+2Lk−Sx01+αL∑k=0pk+2!k!Vak+2Lk=0,(44)Sy0∑k=0pk+1!k!Wak+1Lk−Iy01+αL∑k=0pk+2!k!Uak+2Lk−Ixy01+αL∑k=0pk+2!k!Vak+2Lk=0,(45)Sy02α∑k=0pk+1!k!Wak+1Lk+1+αL∑k=0pk+2!k!Wak+2Lk−Iy03α1+αL∑k=0pk+2!k!Uak+2Lk+1+αL2∑k=0pk+3!k!Uak+3Lk−Ixy03α1+αL∑k=0pk+2!k!Vak+2Lk+1+αL2∑k=0pk+3!k!Vak+3Lk=0,(46)Sx0∑k=0pk+1!k!Wak+1Lk−Ixy01+αL∑k=0pk+2!k!Uak+2Lk−Ix01+αL∑k=0pk+2!k!Vak+2Lk=0,(47)Sx02α∑k=0pk+1!k!Wak+1Lk+1+αL∑k=0pk+2!k!Wak+2Lk−Ixy03α1+αL∑k=0pk+2!k!Uak+2Lk+1+αL2∑k=0pk+3!k!Uak+3Lk−Ix03α1+αL∑k=0pk+2!k!Vak+2Lk+1+αL2∑k=0pk+3!k!Vak+3Lk=0,(48)GJ0∑k=0pk+1!k!ϕak+1Lk−EIω05α1+αL∑k=0pk+2!k!ϕak+2Lk+1+αL2∑k=0pk+3!k!ϕak+3Lk=0,(49)∑k=0pk+2!k!ϕak+2Lk=0.AssumeWa1=c1,Ua2=c2,Ua3=c3,Va2=c4,Va3=c5,ϕa2=c6,ϕa3=c7 and substitute the boundary conditions expressed into the governing equations rewritten in the iterative format using DTM; one of them is presented in equation (41); the transformed function Wak, Uak, Vak, and ϕak can be obtained with iterative process. Consequently, exact solutions of Wa, Ua, Va, and ϕa can be expressed as the sum of all the terms of the corresponding transformed function. The natural frequencies and mode shapes can be determined by applying the boundary conditions and solving a transcendental equation to obtain the unknown ci, i=1,…,7, which is yielded as(50)H11pH12p⋯H17pH21pH22p⋯H27p⋯⋯⋯⋯H71pH72p⋯H77pc1c2⋯c7=0,where Hijpi,j=1,2…7 are polynomials about unknown frequency ωn.Therefore, the eigenvalues can be calculated by taking the determinant of the matrix to be equal to zero:(51)H11pH12p⋯H17pH21pH22p⋯H27p⋯⋯⋯⋯H71pH72p⋯H77p=0.Solving equation (51), the modal frequencies ωn can be obtained. It is necessary to take more terms of Wak, Uak, Vak, and ϕak to evaluate the low order mode frequencies with satisfied precision. This procedure using DFM is different from the CFEM; the corresponding shape function is semianalytical solutions of the governing equations. The iterative calculation for the terms of transformed function Wak, Uak, Vak, and ϕak and modal frequencies of rotor blade can be quickly solved with a symbolic computational software.


## 3.2. Spectral Finite Element Formulation

Similar to the conventional finite element, a dynamic stiffness matrix can be defined to relate the nodal displacements with nodal forces. Figure2 shows the tapered spectral finite element for a rotor blade, in which both nodal displacement and force vectors are defined in the frequency domain, denoted by u^ and F^, respectively(52)u^=U1V1W1ϕ1U2V2W2ϕ2,F^=My1Mx1N1Mt1My2Mx2N2Mt2.Figure 2
Spectral finite element of a tapered thin-walled rotor blade.The nodal force vector can be related with the nodal displacement using a dynamic stiffness matrix,D, which is defined as the following [39]:(53)F^=Dωnu^.Consequently, such standard matrix representation of elemental dynamic stiffness can be used to assemble elements as used in the CFEM. By applying appropriate boundary conditions, we can conduct dynamic analysis of taped thin-walled rotor blade.


## 4. Results and Discussions

In this section, the spectral finite element model is comprehensively validated by comparing modal frequency and mode shape predictions of a tapered cylinder beam and a wind turbine by the ANSYS. Mathematica developed by Wolfram Research, Inc., is used to carry out all calculations in our spectral finite element model.


### 4.1. Tapered Cylinder Beam Analysis

The proposed SFEM is validated by performing modal analysis of a tapered cylinder thin-walled beam. Because of the symmetrical cross section, the warping is not included in its governing equation. Only the flexural modes inx direction are discussed without the coupled axial-bending-torsion effects. Such beam is actually modeled as a classic Euler beam with cantilever boundary conditions. The cylinder thin-walled beam is linearly tapered and its parameters are(54)Iz=I01+αz4,Az=A01+αz2,ρ=2540kg/m3,E=1.7×1010Pa,α=−0.033,μ0=0.6m,L=20m.In the DFM, 50 terms of transform functions are included to compute up to the fourth natural frequency with five-digit precision. The ANSYS software is employed to validate our predictions. Total 16,193 Shell 63 elements were used to obtain comparable accuracy. In our spectral finite element model, one single element is needed. The first four flexural modal frequencies are listed in Table1, SFEM analysis results agree well with ANSYS predictions. The maximum error is less that 1.5%. Figures 3 and 4 show the first four flexural mode shape plots based on the SFEM and simulation results with ANSYS software. Extracting mode displacement results of nodes in ANSYS and drawing the normalized mode shapes, our mode shapes agree well with ANSYS results.Table 1
Comparison of the first four flexural modal frequencies of a tapered cylinder beam.
Order
First mode
Second mode
Third mode
Fourth mode
SFE (Hz)
0.870
3.629
8.631
16.551
ANSYS (Hz)
0.876
3.658
8.716
16.762
Error (%)
0.68
0.79
0.98
1.26Figure 3
First four flexural mode shapes of a tapered cylinder beam using SFEM.Figure 4
The first four flexural mode shapes of a tapered cylinder beam using ANSYS. (a) First mode shape. (b) Second mode shape. (c) Third mode shape. (d) Fourth mode shape.
(a)
(b)
(c)
(d)


### 4.2. Tapered Rotor Blade of Wind Turbine

A tapered rotor blade of a wind turbine with DU40-A17 airfoil shown in Figure1(b) is analyzed with the proposed SFEM in this paper. The blade chord is linearly varied along z axis, as defined in equation (24). The root of the wind turbine blade is fixed to the impeller. Therefore, cantilever boundary conditions are applied, as defined in equations (22) and (23), and the warping is prevented in z=0 and free in z=L [10]. The sectional properties of this single-celled rotor blade calculated using Matlab [38] and material properties are shown below:(55)ρ=2540/m3,E=1.7×1010Pa,G=7.08×109Pa,α=−0.05,L=20m,A0=2.27×10−2m2,J0=4.33×10−3m3,Ix0=1.69×10−3m4,Iy0=1.31×10−2m4,Sx0=−1.00×10−4m3,Sy0=−1.02×10−2m3,Ixy0=−4.42×10−5m4,Iω0=8.65×10−6m6.Total 30,000 Shell 63 elements are used to model the rotor blade, as shown in Figure5. One single spectral finite element is needed. Table 2 shows the modal frequency results. Our results agree with ANSYS predictions. The maximum error is less than 9%. It is noted that such errors are larger than the errors of the first four flexural modes of the above tapered cylinder beam. Such extra error is caused by the coupled axial-bending-torsion effect in the tapered rotor blade.Figure 5
FE model in ANSYS of the tapered wind turbine blade.Table 2
Comparison of modal frequencies of the tapered wind turbine blade.
Order
First mode
Second mode
Third mode
Fourth mode
Mode shape
Flapwise
Edgewise
Flapwise
Flapwise
SFE (Hz)
1.060
2.392
3.960
5.131
ANSYS (Hz)
1.016
2.240
4.299
5.559
Error (%)
4.32
6.78
8.56
8.34The coupled mode shapes are shown in Figure6, where “Rz” denotes the torsional displacement. The flapwise bending mode is dominated in the first, the third, and the fourth mode. The edgewise bending mode is dominated in the second mode. Due to the coupled motion, both axial and torsion modes appear in these bending dominated modes. It is crucial to include all motions in the dynamic analysis of a rotor blade. The first four mode shapes of the wind turbine blade with ANSYS software are shown in Figure 7, but the coupling among axial, bending, and torsion motion is complicated to be drawn like the mode shape plots in Figure 6. It is noted that the dominated flapwise and edgewise mode shapes agree well with our mode shapes.Figure 6
The first four mode shapes of the tapered wind turbine blade using SFEM. (a) First mode shape. (b) Second mode shape. (c) Third mode shape. (d) Fourth mode shape.
(a)
(b)
(c)
(d)Figure 7
The first four mode shapes of the tapered wind turbine blade using ANSYS. (a) First mode shape. (b) Second mode shape. (c) Third mode shape. (d) Fourth mode shape.
(a)
(b)
(c)
(d)


## 4.1. Tapered Cylinder Beam Analysis

The proposed SFEM is validated by performing modal analysis of a tapered cylinder thin-walled beam. Because of the symmetrical cross section, the warping is not included in its governing equation. Only the flexural modes inx direction are discussed without the coupled axial-bending-torsion effects. Such beam is actually modeled as a classic Euler beam with cantilever boundary conditions. The cylinder thin-walled beam is linearly tapered and its parameters are(54)Iz=I01+αz4,Az=A01+αz2,ρ=2540kg/m3,E=1.7×1010Pa,α=−0.033,μ0=0.6m,L=20m.In the DFM, 50 terms of transform functions are included to compute up to the fourth natural frequency with five-digit precision. The ANSYS software is employed to validate our predictions. Total 16,193 Shell 63 elements were used to obtain comparable accuracy. In our spectral finite element model, one single element is needed. The first four flexural modal frequencies are listed in Table1, SFEM analysis results agree well with ANSYS predictions. The maximum error is less that 1.5%. Figures 3 and 4 show the first four flexural mode shape plots based on the SFEM and simulation results with ANSYS software. Extracting mode displacement results of nodes in ANSYS and drawing the normalized mode shapes, our mode shapes agree well with ANSYS results.Table 1
Comparison of the first four flexural modal frequencies of a tapered cylinder beam.
Order
First mode
Second mode
Third mode
Fourth mode
SFE (Hz)
0.870
3.629
8.631
16.551
ANSYS (Hz)
0.876
3.658
8.716
16.762
Error (%)
0.68
0.79
0.98
1.26Figure 3
First four flexural mode shapes of a tapered cylinder beam using SFEM.Figure 4
The first four flexural mode shapes of a tapered cylinder beam using ANSYS. (a) First mode shape. (b) Second mode shape. (c) Third mode shape. (d) Fourth mode shape.
(a)
(b)
(c)
(d)


## 4.2. Tapered Rotor Blade of Wind Turbine

A tapered rotor blade of a wind turbine with DU40-A17 airfoil shown in Figure1(b) is analyzed with the proposed SFEM in this paper. The blade chord is linearly varied along z axis, as defined in equation (24). The root of the wind turbine blade is fixed to the impeller. Therefore, cantilever boundary conditions are applied, as defined in equations (22) and (23), and the warping is prevented in z=0 and free in z=L [10]. The sectional properties of this single-celled rotor blade calculated using Matlab [38] and material properties are shown below:(55)ρ=2540/m3,E=1.7×1010Pa,G=7.08×109Pa,α=−0.05,L=20m,A0=2.27×10−2m2,J0=4.33×10−3m3,Ix0=1.69×10−3m4,Iy0=1.31×10−2m4,Sx0=−1.00×10−4m3,Sy0=−1.02×10−2m3,Ixy0=−4.42×10−5m4,Iω0=8.65×10−6m6.Total 30,000 Shell 63 elements are used to model the rotor blade, as shown in Figure5. One single spectral finite element is needed. Table 2 shows the modal frequency results. Our results agree with ANSYS predictions. The maximum error is less than 9%. It is noted that such errors are larger than the errors of the first four flexural modes of the above tapered cylinder beam. Such extra error is caused by the coupled axial-bending-torsion effect in the tapered rotor blade.Figure 5
FE model in ANSYS of the tapered wind turbine blade.Table 2
Comparison of modal frequencies of the tapered wind turbine blade.
Order
First mode
Second mode
Third mode
Fourth mode
Mode shape
Flapwise
Edgewise
Flapwise
Flapwise
SFE (Hz)
1.060
2.392
3.960
5.131
ANSYS (Hz)
1.016
2.240
4.299
5.559
Error (%)
4.32
6.78
8.56
8.34The coupled mode shapes are shown in Figure6, where “Rz” denotes the torsional displacement. The flapwise bending mode is dominated in the first, the third, and the fourth mode. The edgewise bending mode is dominated in the second mode. Due to the coupled motion, both axial and torsion modes appear in these bending dominated modes. It is crucial to include all motions in the dynamic analysis of a rotor blade. The first four mode shapes of the wind turbine blade with ANSYS software are shown in Figure 7, but the coupling among axial, bending, and torsion motion is complicated to be drawn like the mode shape plots in Figure 6. It is noted that the dominated flapwise and edgewise mode shapes agree well with our mode shapes.Figure 6
The first four mode shapes of the tapered wind turbine blade using SFEM. (a) First mode shape. (b) Second mode shape. (c) Third mode shape. (d) Fourth mode shape.
(a)
(b)
(c)
(d)Figure 7
The first four mode shapes of the tapered wind turbine blade using ANSYS. (a) First mode shape. (b) Second mode shape. (c) Third mode shape. (d) Fourth mode shape.
(a)
(b)
(c)
(d)


## 5. Conclusions

In this paper, a spectral finite element model is developed to efficiently and accurately predict the dynamic behavior of tapered thin-walled rotor blade. We first deduced the fully coupled governing equations using Hamilton’s principle, in which axial, bending, and torsion motions are included. Subsequently, the DTM was applied to obtain the semianalytical solutions in order to formulate the spectral finite element. Modal analyses of a tapered cylinder beam and a tapered thin-walled wind turbine rotor blade were conducted to validate the spectral finite element model. Modal frequency results agreed well with the ANSYS predictions; the maximum error is less than 9% for the example wind turbine blade. Only one single spectral finite element is needed because the interpolation functions are duplicated from the exact semianalytical solutions. The coupled modal characteristics among flapwise and edgewise bending, axial, and torsion motion were captured. These modal frequency and mode shape results are beneficial to carry out rotor blade design and performance analysis for wind turbine and helicopter.



---

*Source: 2174209-2019-01-02.xml*