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

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. )en, the differential transformmethod (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. )ese modal frequency andmode shape results are important to carry out design and performance evaluation of the tapered thin-walled structures.


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. erefore, 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. e need of onedimensional (1D) beam finite element model for tapered thin-walled beam still remains as an attractive research topic. e 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 closedsection 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 flexuraltorsional 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 thinwalled structure, including axial, bending, torsion, and torsional-related warping, as well as the anisotropic thinwalled beam modeling of aircraft wing. Recently, improved 1D beam models are developed for thin-walled beam by Vo and Lee [10][11][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][16][17][18][19][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. e 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]. e SFEM is developed in the frequency domain [27,28]. e 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. e 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]. e 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. en, 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. e paper is organized as follows. Section 2 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.

Modeling of Tapered Thin-Walled
Rotor Blades

Kinematics.
A generic single-cell tapered thin-walled rotor blade is shown in Figure 1, in which the airfoil shape is DU40-A17 of a NREL 5 MW wind turbine [36]. e 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). e first set is a global Cartesian coordinate system oxyz. e z axis is the axial direction of the rotor blade, and its position is placed through the shear centers of the cross sections herein. e x and y axes are the transverse coordinates of the closed cross section. e 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. erefore, 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) e contour does not deform in its own plane (2) e 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) e Kirchhoff-Love assumption is applied to the thin wall Assumption 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]: 2 Shock and Vibration As shown in Figure 1(b), the associated parameters r(s) and q(s) are given by where θ(s) is the angle between the x axis and the local s axis at the point s on the perimeter. e corresponding out-of-plane shell displacement w can be obtained based on assumption 2. Neglecting the transverse shearing effects on the deformation of the blade, the shear strain c sz of the middle surface can be expressed as 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: where dx � ds sin θ(s) and dy � ds cos θ(s).
According to the Kirchhoff-Love shell theory in assumption 3, the local displacements at any point P on the shell wall can be expressed as e strains in the shell wall are given by Finally, the strains in the shell wall can be obtained by substituting equation (8) into equation (7) as shown below:

Governing
Equations. e total strain energy K is given by Substituting equation (9) into equation (10), the variation of strain energy can be rewritten in terms of general force and displacement components where N z , M y , M x , M t , and M ω are the axial force, bending moments, torsional moment, and warping bimoment, respectively. ey are defined by integrating over the cross section as e final expression of these force components are given below: Corresponding sectional properties are varied along z axis for a tapered thin-walled beam. ese parameters at an arbitrary cross section can be simplified for a thin-walled beam and are shown below: 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: ere 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. e other parameters are calculated by using Shi's finite segment method [38]. e kinematic energy is given by e variation of the kinematic energy is shown below: Applying Hamilton's principle, the governing equations and associated boundary conditions are obtained. e coupled axial-bending-torsion equations are e blade is a cantilever beam, corresponding boundary conditions are given by z � 0 : W(0) � 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 4 Shock and Vibration where α is the taper ratio, α � (μ L /μ 0 − 1)/L and μ L and μ 0 are the chord length of its two ends, respectively. e sectional properties at any cross section can be rewritten as follows in order to account for the taper: Finally, the governing equations (18)-(21) can be modified to account for the variation of sectional properties:

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. en, 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.

Semianalytical Solutions.
Assuming harmonic responses, the displacements U, V, W, and ϕ in the time domain can be written in the following forms [39]: where U a , V a , W a , and ϕ a are the steady state frequency dependent displacement solutions. ω n is the angular excitation frequency. e governing equations (26)- (29) can be rewritten as Shock and Vibration ρS x0 ω 2 n U a − ρS y0 ω 2 n V a − ρω 2 n I x0 + I y0 (1 + αz)ϕ a + ρI ω0 ω 2 n 5α(1 + αz) 2 Note that W ′″ a can be obtained by calculating the differential of W ″ a from equation (31). It yields Substituting equation (35) into equations (32) and (33), the corresponding governing equations can be rewritten as

Definition in the DFM.
e 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. e basic definitions and the application procedure of DFM are presented in the following.
In DFM, the differential transform of the function f(z) is defined as where f(z) is the original function and F(k) is the transformed function. e inverse transformation can be expressed as where p denotes the finite terms F(k) of the transformed function, which depends on the convergence rate of the modal frequency prediction of rotor blade. eorems that are used in the transformation of the governing differential equations and the boundary conditions of rotor blade are expressed as [34]

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: 6 Shock and Vibration Furthermore, the boundary conditions in equation (22) can be rewritten as Considering equation (13), the boundary conditions shown in equation (24) are yielded as the following equations: Solving equation (51), the modal frequencies ω n can be obtained. It is necessary to take more terms of W a (k), U a (k), V a (k), and ϕ a (k) to evaluate the low order mode frequencies with satisfied precision. is procedure using DFM is different from the CFEM; the corresponding shape function is semianalytical solutions of the governing equations. e iterative calculation for the terms of transformed function W a (k), U a (k), V a (k), and ϕ a (k) and modal frequencies of rotor blade can be quickly solved with a symbolic computational software.

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. Figure 2 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 e nodal force vector can be related with the nodal displacement using a dynamic stiffness matrix, D, which is defined as the following [39]: 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.

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.

Tapered Cylinder Beam Analysis.
e 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 in x direction are discussed without the coupled axial-bending-torsion effects. Such beam is actually modeled as a classic Euler beam with cantilever boundary conditions. e cylinder thin-walled beam is linearly tapered and its parameters are In the DFM, 50 terms of transform functions are included to compute up to the fourth natural frequency with five-digit precision. e 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. e first four flexural modal frequencies are listed in Table 1, SFEM analysis results agree well with ANSYS predictions. e 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.

Tapered Rotor Blade of Wind Turbine.
A tapered rotor blade of a wind turbine with DU40-A17 airfoil shown in Figure 1(b) is analyzed with the proposed SFEM in this paper. e blade chord is linearly varied along z axis, as defined in equation (24). e root of the wind turbine blade is fixed to the impeller. erefore, 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]. e sectional properties of this single-celled rotor blade calculated using Matlab [38] and material properties are shown below: (55) Total 30,000 Shell 63 elements are used to model the rotor blade, as shown in Figure 5. One single spectral finite element is needed. Table 2 shows the modal frequency results. Our results agree with ANSYS predictions. e 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. e coupled mode shapes are shown in Figure 6, where "Rz" denotes the torsional displacement. e flapwise bending mode is dominated in the first, the third, and the fourth mode. e 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. e 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.

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. e coupled modal characteristics among flapwise and edgewise bending, axial, and torsion motion were captured. ese modal frequency and mode shape results are beneficial to carry out rotor blade design and performance analysis for wind turbine and helicopter.
Data Availability e airfoil data of DU40 and FE model of the tapered blade, as well as the program to calculate the cross section properties and the eigenvalue, are published in the figshare database (https://figshare.com/articles/Airfoil_data_cross-section_parameter_calculation_and_DTM-SFEM_program_ and_FE_model_of_the_tapered_blade/7001243).

Conflicts of Interest
e authors declare that they have no conflicts of interest.