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.

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 [1012], 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 [1520]. 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. 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.

2. Modeling of Tapered Thin-Walled Rotor Blades

2.1. 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]. 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).

The first set is a global Cartesian coordinate system . The axis is the axial direction of the rotor blade, and its position is placed through the shear centers of the cross sections herein. The and axes are the transverse coordinates of the closed cross section. The second set is a local shell wall coordinate system , where is a contour coordinate measured along the middle surface of the shell wall and 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 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 and at any point in the local shell wall coordinate system can be expressed in terms of displacements , in the and directions, respectively, and the rotation angle about axis in the global coordinate system [11, 15]:

As shown in Figure 1(b), the associated parameters and are given bywhere 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 displacement can be obtained based on assumption 2. Neglecting the transverse shearing effects on the deformation of the blade, the shear strain 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 , the axial displacement can be determined as shown below:where and .

According to the Kirchhoff–Love shell theory in assumption 3, the local displacements at any point on the shell wall can be expressed as

The 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:

2.2. Governing Equations

The total strain energy is given by

Substituting equation (9) into equation (10), the variation of strain energy can be rewritten in terms of general force and displacement componentswhere are the axial force, bending moments, torsional moment, and warping bimoment, respectively. They are defined by integrating over the cross section as

The final expression of these force components are given below:

Corresponding sectional properties are varied along 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:

It has been discussed by Gjelsvik [4] that the following parameters are zero when the pole of the local sectional coordinate system is identical to the shear center:

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

The variation of the kinematic energy is shown below:

Applying Hamilton’s principle, the governing equations and associated boundary conditions are obtained. The coupled axial-bending-torsion equations are

The blade is a cantilever beam, corresponding boundary conditions are given by

In this study, it is assumed that the chord length of rotor blade is uniformly tapered along axis. It is defined bywhere is the taper ratio, and and 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:

Finally, the governing equations (18)–(21) can be modified to account for the variation of sectional properties:

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 displacements in the time domain can be written in the following forms [39]:where are the steady state frequency dependent displacement solutions. is the angular excitation frequency.

The governing equations (26)–(29) can be rewritten as

Note that can be obtained by calculating the differential of from equation (31). It yields

Substituting equation (35) into equations (32) and (33), the corresponding governing equations can be rewritten aswhere .

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 function is defined aswhere is the original function and is the transformed function.

The inverse transformation can be expressed aswhere denotes the finite terms 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]

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:

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:

Assume 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 , , , and can be obtained with iterative process. Consequently, exact solutions of , , , and 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 , , which is yielded aswhere are polynomials about unknown frequency .

Therefore, the eigenvalues can be calculated by taking the determinant of the matrix to be equal to zero:

Solving equation (51), the modal frequencies can be obtained. It is necessary to take more terms of , , , and 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 , , , and 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. 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 and , respectively

The nodal force vector can be related with the nodal displacement using a dynamic stiffness matrix, , 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.

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 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. The 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. 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 Table 1, 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.

4.2. 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. The blade chord is linearly varied along 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 and free in [10]. The sectional properties of this single-celled rotor blade calculated using Matlab [38] and material properties are shown below:

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

The coupled mode shapes are shown in Figure 6, 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.

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.

Data Availability

The 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

The authors declare that they have no conflicts of interest.

Acknowledgments

This study was supported by the National Natural Science Foundation of China (nos. 11572125, 51575176, and 51775182).