Abstract
Based on free interface component modal synthesis method, the free vibration behavior of wind turbine structures is investigated. The wind turbine structure is divided into three parts including tower, wheel hubcabin and rotor. The tower is modeled as an isotropic metal cantilever beam, the blade as thinwalled composite beam and the wheel hubcabin as a rigid body due to its large extensional stiffness, bending stiffness and torsion stiffness compared with tower and blades. The displacements of the blades are described by thinwalled composite beam theory. Galerkin’s method is used to discretize blades and tower. Employing Lagrange method, the motion equations of blades are derived and then stiffness and mass matrices are obtained. The natural frequencies and mode shapes of the wind turbine structure are predicted by numerical simulations. Numerical results using the present model are validated by ANSYS software results.
1. Introduction
Large horizontal axis wind turbines are sophisticated dynamic systems in which structural dynamic and aeroelastic stability problems can be frequently suffered from. These problems appear because of the highly variable turbulent atmospheric wind environment that flexible composite blade system encounters, as it rotates. In order to further improve the dynamical performance of wind power generators, a comprehensive model is therefore essential.
Currently, the most widely used structural modeling method of wind turbine is multibody dynamical method. Several multibody dynamical formulations have been presented for the analysis of horizontal axis wind turbines [17]. Typically, flexible components, such as tower and rotor blades, are discretized by an FE method. In general, in order to provide an accurate system model, many degrees of freedom for each component need to be involved. Thus, the total degrees of freedom will increase significantly, which will in turn increase the computational effort, especially for such timedependent systems as large wind power generator systems. Substructure or component modal synthesis is useful to model large and complex mechanical systems where the exact solution is not available or numerical techniques are inefficient. However, to the best of the author’s knowledge, no literature on structural dynamic modelling of a wind turbine by component mode synthesis is cited, though a recent study has been found in the dynamical analysis of the wind turbines FE based simulations for rotor blades [6].
The main problem of in the design of wind turbine is load carrying capacity and the fatigue resistance due to dynamic instabilities. As first part of the dynamical design of wind turbine, the structural modelling (no aerodynamic loads is included) and analysis of the basic dynamic properties, such as natural frequencies and mode shapes are of vital importance for further investigation of stability problems and structural integrity of the entire wind turbine. The goal of this paper is to present a simpler and less computationally expensive model in order to investigate dynamic characteristics of wind turbine. Free interface component modal synthesis method is used to calculate the natural frequency of wind turbine. The tower is modeled as an isotropic metal cantilever beam, the rotor blades as composite thinwalled closed section beam. A composite thinwalled beam theory, which is an asymptotically correct theory referred to as variational asymptotically method (VAM) [8] is applied to obtain the displacement field in the rotor. Galerkin’s approach is used for the discretization in space for blades and tower. By using the continuity conditions of the deformation, mass and stiffness matrices of these flexible components are assembled to get mass and stiffness matrices of whole wind turbine system. Natural frequencies of wind turbine are evaluated from the free vibration solution. The finite element analysis software, ANSYS is also used to calculate the natural frequencies of wind turbine in order to examine validity of the present model.
2. Wind turbine modeling
The deformations of cabin and wheel hub are very smaller than the blade and tower. So it’s reasonable to model the cabin and wheel hub as rigid body, whereas blade and tower as flexible body. The structural model and coordinate systems are shown in Fig. 1. Inertia coordinate system $\left(OXYZ,\mathrm{}\overrightarrow{I},\overrightarrow{J},\mathrm{}\overrightarrow{K}\right)$, which is fixed, and coordinate origin is located at the right center of wheel hub before deformation. Tower coordinate system $\left({O}_{t}{x}_{t}{y}_{t}{z}_{t},\mathrm{}\overrightarrow{i},\mathrm{}\overrightarrow{j},\mathrm{}\overrightarrow{k}\right)$, which describes the tower motion, and coordinate origin is located at the root of tower. And it doesn’t move during vibration. Cabin coordinate system $\left({O}_{j}{x}_{j}{y}_{j}{z}_{j},\overrightarrow{i},\mathrm{}\overrightarrow{j},\mathrm{}\overrightarrow{k}\right)$, which describes cabin motion. The coordinate origin is located at the center of cabin, and just at the top of tower. It moves along with the tower wheel hub coordinate system $\left({O}_{s}{x}_{s}{y}_{s}{z}_{s},{\overrightarrow{i}}_{s},\mathrm{}{\overrightarrow{j}}_{s},\mathrm{}{\overrightarrow{k}}_{s}\right)$, which describes the motion of wheel hub. The coordinate origin is located at the wheel hub center. And it moves with the cabin and tower. Blade coordinate system $\left({O}_{s}xyz,\overrightarrow{i},\mathrm{}\overrightarrow{j},\mathrm{}\overrightarrow{k}\right)$, which describes the blade motion. Its coordinate origin is located at the wheel hub and moved with the cabin and tower.
Based on the structural characteristics of wind turbine, it is divided into three parts: rotor, wheel hub and cabin, tower. This is convenient to apply modal synthesis method.
Fig. 1Wind turbine structural model and coordinate systems
Fig. 2Inertia coordinate system and wheel hub coordinate system after deformation
3. Free vibration characteristics of subsystems
3.1. Free vibration characteristics of rotor
The wheel hub coordinate system and blade coordinate system are moveable because of the displacements of tower and cabin as shown in Fig. 2. ${q}_{x}$, ${q}_{y}$, ${q}_{z}$ are displacements and ${\phi}_{x}$, ${\phi}_{y}$, ${\phi}_{z}$ twist angles of wheel hub in inertia coordinate system.
Besides all above rigid displacements, the blades have elastic deformations. Let ${U}_{1i}$, ${U}_{2i}$ and ${U}_{3i}$ represent the average displacement along $x$, $y$, $z$ directions while ${\phi}_{i}$ represents the twist angle about $x$axis, and the subscript $i$ ($i=$ 1, 2, 3) identifies the ith blade. So the rotor has 6 rigid displacements and 12 elastic displacements.
3.1.1. Transformation matrix from inertia coordinate system to blade coordinate system
In order to derive the kinetic energy of rotor, the rotate transformation matrix from inertia coordinates to blade coordinates must be used.
The inertia coordinate rotates ${\phi}_{y}$ about $y$axis, ${\phi}_{x}$ about $x$axis and ${\phi}_{z}$ about $z$axis to the coordinate system of first blade. It rotates same angle about $x$, $y$axis to the coordinate system of second and the third blades. While rotate $z+\text{120\xb0}$ and $z120\xb0$ about $z$axis to the coordinate system of second and third blades, respectively.
Letting $\mathrm{sin}{\phi}_{x}\cong {\phi}_{x}$, $\mathrm{sin}{\phi}_{y}\cong {\mathrm{\phi}}_{\mathrm{y}}$, $\mathrm{sin}{\phi}_{z}\cong {\phi}_{z}$, the transformation matrix from inertia to the blades’s coordinate system is:
3.1.2. Kinetic energy density of rotor
The displacement vector of any point on the first blade can be expressed in terms of the inertia coordinate system:
where ${u}_{11}$, ${u}_{21}$, ${u}_{31}$ are the displacements of any point on the crosssection of the first blade in $x$, $y$, $z$ direction. Based on Ref. [8], the displacement fields are assumed:
${u}_{21}={U}_{21}z{\phi}_{1},$
${u}_{31}={U}_{31}+y{\phi}_{1}.$
Based on the displacement vector and coordinate transformation, the velocity of the first blade is expressed as:
The kinetic energy of per unit length is:
By combining the transformation matrix with Eq. (4), Eq. (6) and Eq. (7), the kinetic energy of the first blade can be obtained. And similarly, the kinetic energy of the second blade ${k}_{2}$ and the third blade ${k}_{3}$ can also be obtained.
Thus the kinetic density of rotor takes the form:
3.1.3. Strain energy density of rotor
Based on Ref. [9], the per unit length strain energy of the first blade takes the form:
where, $\left\{\delta \right\}=\left\{{U}_{11}^{\text{'}}{\phi}_{1}^{\text{'}}{U}_{31}^{\text{'}\text{'}}{U}_{21}^{\text{'}\text{'}}\right\}$.
In Eq. (9), the first term is the strain energy due to deformation of the composite blade, the second and third terms are the strain energy due to rotating motion of the blade, where, ${\left[C\right]}_{4\times 4}$ is a 4×4 symmetric stiffness matrix, ${m}_{c}=\oint \rho h\left(s\right)ds$ is the mass of per unit length, $h\left(s\right)$ is the thickness of the blade at any point, ${\left(\mathrm{}\mathrm{}\right)}^{\mathrm{\text{'}}}$ and ${\left(\mathrm{}\mathrm{}\right)}^{\text{'}\text{'}}$ denote first and second differentiation with respect to $x$, respectively, $\left(\dot{}\right)$ and $\left(\ddot{}\right)$ denote first and second differentiation with respect to $t$, respectively.
The strain energy density ${U}_{2}$ and ${U}_{3}$ of second and third blades can also be obtained just by substituting the displacement variables into Eq. (5) with their corresponding variables.
So the strain energy density of the rotor can be written as:
3.1.4. Gravitational potential energy density of rotor
The gravity acts on the blades along the inertia negative $X$ axis. The gravitational potential energy density of the first blade takes the following form:
where, ${X}_{1}$ can be found in Eq. (4).
Next, one can obtain the gravitational potential energy density of the second blade ${V}_{g2}$ and the third blade ${V}_{g3}$.
Thus gravitational potential energy density of the rotor is:
The Lagrange density of rotor [10] can be written as:
3.1.5. Mass matrix and stiffness matrix of rotor
The mode shape of rigid motion variable ${q}_{x}$, ${q}_{v}$, ${q}_{z}$ and ${\phi}_{x}$, ${\phi}_{v}$, ${\phi}_{z}$ is assumed to be $\overrightarrow{E}=(\mathrm{1,1},\dots ,1)$. And elastic deformation variable is assumed to take the standard nonrotating, noncoupling, uniform cantilever beam mode shape functions as:
${\theta}_{j}=\sqrt{2}\mathrm{sin}\left(\frac{{\gamma}_{j}x}{L}\right),\mathrm{cos}\left({\beta}_{j}\right)\mathrm{cos}h\left({\beta}_{j}\right)=1,$
${\lambda}_{j}=\frac{\left(\mathrm{cos}{\beta}_{j}+\mathrm{cosh}{\beta}_{j}\right)}{\left(\mathrm{sin}{\beta}_{j}+\mathrm{sinh}{\beta}_{j}\right)},$
${\psi}_{j}=\mathrm{cos}\left(\frac{{\beta}_{j}x}{L}\right)\mathrm{cosh}\left(\frac{{\beta}_{j}x}{L}\right)+{\lambda}_{j}\left(\mathrm{sin}\left(\frac{{\beta}_{j}x}{L}\right)\mathrm{sinh}\left(\frac{{\beta}_{j}x}{L}\right)\right),j=1,2,\dots ,N,$
where, $N$ is number of mode shapes.
Substitution the mode shape function Eq. (14) into Eq. (13), then integrating along the blade length, the Lagrange function can be obtained. Based on the Lagrange equation, the elements of mass matrix ${M}_{y}$ and stiffness matrix ${K}_{y}$ of the rotor are just the coefficients of corresponding variables, so the ${M}_{y}$ and ${K}_{y}$ can be derived directly according to the Lagrange function.
3.2. Free vibration characteristics of wheel hubcabin and tower
The tower, wheel hubcabin are considered as simple isotropic beams. So their vibrations dispalcements can be discribed based on the classical extension, bending and torsion beam theory. And the tower is a cantilever beam, and wheel hubcabin can be simplified as a freefree beam.
The tower has 4 degrees of freedom by ${w}_{tx}$, ${w}_{tv}$, ${w}_{tz}$, ${\phi}_{tx}$, where, ${w}_{tx}$ denotes displacements along ${x}_{t}$, ${w}_{tv}$ denotes deflection along ${y}_{t}$, ${w}_{tz}$ denotes deflection along ${z}_{t}$ and ${\phi}_{tx}$ denotes torsion about ${x}_{t}$. And its mode shape is same as in Eq. (14). While the wheel hub and cabin have 4 elastic displacements and 6 rigid displacements caused by the deformation of tower.
The motions of cabin and wheel in every direction can be expressed as:
Displacement along ${x}_{j}$ direction:
Displacement along ${z}_{j}$ direction:
Displacement along ${y}_{j}$ direction:
Rotation about $x$:
where, the mode shape of rigid motion is also supposed to be $\overrightarrow{E}=\left(\mathrm{1,1},\dots ,1\right)$, and the mode shapes of freefree beam can be represented in the form:
${\theta}_{ji}=\sqrt{2}cos\left(\frac{i\pi}{{L}_{j}}x\right),{\gamma}_{i}=\frac{\mathit{cosh}{\beta}_{i}\mathit{cos}{\beta}_{i}}{\mathit{sin}h{\beta}_{i}\mathit{sin}{\beta}_{i}},1\mathit{cos}{\beta}_{i}\mathit{cosh}{\beta}_{i}=0,i=\mathrm{1,2},\dots ,N,$
where, $N$ is number of mode shape.
Based on Galerkin’s method, one can obtain the mass matrix and stiffness matrix of tower, wheel hub and cabin, which are denoted by ${M}_{j}$, ${K}_{j}$, ${M}_{t}$, ${K}_{t}$.
4. Free vibration characteristic of wind turbine
By assembling the subsystem’s mass and stiffness matrices, the whole wind turbine mass and stiffness matrices take the form:
where, ${M}^{\text{'}}$ and ${K}^{\text{'}}$ are 32$N$×32$N$ symmetric matrices. The 32 variables are not independent. The deformations are continuous at the interface of tower and cabin, and at the interface of wheel hub and rotor. So the following relations have to be fulfilled:
${\phi}_{tx}\left(H,t\right)={W}_{jv}^{\text{'}}\left(0,t\right),{W}_{tz}^{\text{'}}\left(H\right)={W}_{jx}^{\text{'}}\left(0\right),{W}_{tv}^{\text{'}}\left(H\right)={\phi}_{jz}\left(0\right),$
${q}_{x}\left(t\right)={W}_{jx}\left({L}_{j},t\right),{q}_{j}\left(t\right)={W}_{jy}\left({L}_{j},t\right),{q}_{z}\left(t\right)={W}_{jz}\left({L}_{j},t\right),$
${\phi}_{x}\left(t\right)={W}_{jy}^{\text{'}}\left({L}_{j},t\right),{\phi}_{y}\left(t\right)={W}_{jx}^{\text{'}}\left({L}_{j},t\right),{\phi}_{z}\left(t\right)={\phi}_{jz}\left({L}_{j},t\right).$
By using Eqs. (21), the transformation matrix $D$ can be formed from which one can reduce ${M}^{\text{'}}$ and ${K}^{\text{'}}$ from 32$N$ variables to 20$N$ variables. Thus, the mass and stiffness matrices of wind turbine described as Eqs. (20) become:
The expressions for ${M}^{\text{'}}$, ${K}^{\text{'}}$, $M$, $K$, and $D$ are not shown for the sake of simplicity. From $M$ and $K$, the generalized eigenvalue problem can be yielded. So the free vibration of wind turbine system can be studied by the eigenvalue equations.
5. Numerical examples
The blades is a thinwalled box beam made of T300/5208 graphite/epoxy composite material, while tower and wheel hubcabin are isotropic hollow cylinder beams. The geometric and material properties of the blades, the tower and the wheel hubcabin are given in Table 1 and Table 2, respectively. These property parameters are chosen from 2.5 MW wind turbine data. Applying Eqs. (22) to obtain the mass and stiffness matrices of wind turbine, the free vibration frequencies of the wind turbine can be calculated.
Table 1The geometric and material properties of the blades
Length $=$ 44.2 m  Ply thickness $=$ 0.05 cm 
Width $=$ 2 m  Depth $=$ 0.4 m 
Number of piles $=$ 100  $\rho =\text{1600}$ kg/m^{3} 
${E}_{11}=$142 GPa  ${E}_{22}={E}_{33}=$9.8 GPa 
${G}_{12}={G}_{13}=$6 GPa  ${G}_{23}=$4.83 GPa 
${\mu}_{12}={\mu}_{13}=\text{0.42}$  ${\mu}_{23}=\text{0.5}$ 
The first five natural frequencies of the wind turbine obtained by the present model in the still state are shown in Table 4. They are compared with the corresponding results of ANSYS software. For ANSYS, SHELL 99 element is chosen to discrete the blade, while BEAM44 is used to discrete the tower, wheel hub and cabin. The connects between the wheel hub and blades are assumed as rigid constraints. Gravity and rotating speed are loaded on the blades. A perfect agreement of numerical results using the present model with ANSYS results can be seen from Table 4.
Table 2The geometric and material properties of the tower
Height $=$ 80 m  Outer diameter $=$ 2.515 m 
Inner diameter $=$ 2.485 m  ${E}_{t}=$200 GPa 
${\mu}_{t}=$0.3  ${\rho}_{\tau}=\text{7800}$ kg/m^{3} 
Table 3The geometric and material properties of the cabinwheel hub
Length $=$ 3 m  Outer diameter $=$ 2.3 m 
Inner diameter $=$ 2 m  ${E}_{j}=$20000 GPa 
${\mu}_{j}$$=$ 0.3  ${\rho}_{j}=\text{7800}$ kg/m^{3} 
Table 4First five natural frequencies (Hz) of the wind turbine with Ω= 0 r/min
Method  1st^{}  2nd  3rd  4th  5th 
Modal Synthesis  0.1271  0.1396  0.1435  0.1689  0.2474 
ANSYS  0.1242  0.1272  0.1426  0.1686  0.2355 
The variation of first five natural frequencies vs. rotating speed of rotor blades is shown in Fig. 3 from the present model with the ones from ANSYS software. As it appears from these figures, natural frequencies increase as the rotating speed is increased, because of the existence of a larger centrifugal forces at a higher rev/min, which is called as the dynamic stiffening. However, it can be noted that the effect of the rotating speed appears to be not significant for the lowermode ones. It can also be observed that the results from the present model in general, agree well, as compared with ANSYS, at lower rotating speed. Fig. 3 shows close agreement between the present model and ANSYS results for the first mode, second mode and fouth mode. And the third mode and fifth mode depart slightly in magnitude at higher rotating speed but still have the same trend.
Fig. 3Variation of first five natural frequencies with rotating speed
The first five mode shapes are shown in Fig. 4 to Fig. 8. The most important dynamic properties of the complete wind turbine structural systems can be seen from these figures, they are: the flap and lag vibrations of rotor blades (see Figs. 45), the vibration interactions of the rotor blades and tower (see Figs. 68).
Fig. 4The first mode shape
Fig. 5The second mode shape
Fig. 6The third mode shape
Fig. 7The fourth mode shape
Fig. 8The fifth mode shape
6. Conclusions
This paper deals with the structural dynamical modeling of wind turbine. An analytical model capable of predicting of the natural frequencies of wind turbine has been proposed based on free interface component modal synthesis method. The present model employs thinwalled composite beam theory that can be used effectively for the analysis of tubular thinwalled composite blades. Flexible subsystems as blades and tower are discretized using Galerkin’s method that is computationally efficient compared to the traditional finite element method when applied to compute the dynamic behavior of wind turbine. Results obtained by the present model are agreement with the ANSYS results. Our study exhibits the capability of component modal synthesis method to model the wind turbine. This work is a part of a project on structural dynamics and aeroelastics of wind turbine, and the model will be extended to include aerodynamic loads arising from the imposed wind field and be used for analysis of stability problems of wind turbine.
References

Stol K. Dynamics modeling and periodic control of horizontalaxis wind turbines. Ph. D. Thesis, University of Colorado at Boulder, Boulder CO, USA, 2001.

Zhao X. Y., Peter M., Wu J. Y. A new multibodymodeling methodology for wind turbine structures using a cardanic joint beam element. Renewable Energy, Vol. 32, Issue 3, 2007, p. 532546.

Lee D. H., Dewey H., Hodges P., Mayuresh J. P. Multiflexiblebody dynamic analysisof horizontal axis wind turbines. Wind Energy, Vol. 5, Issue 4, 2002, p. 281300.

MaiXer P. A differentialgeometric approach to the multibody system dynamics. Journal of Applied Mathematics and Mechanics, Vol. 71, Issue 4, 1991, p. 116119, (in German).

Schiehlen W. Advanced Multibody System DynamicsSimulation and Software Tools. London, KluwerAcademic Publishers, 1993.

HolmJorgensen K.,Nielsen S. R. K. A component mode synthesis algorithm for multibody dynamics ofwind turbines. Journal of Sound and Vibration, Vol. 326, Issue 35, 2009, p. 53767.

Wang J. H., Qin D.T., Lim T. C. Dynamic analysis of horizontal axis wind turbine by thinwalled beam theory. Journal of Sound and Vibration, Vol. 329, Issue 17, 2010, p. 35653586.

Berdichevsky V. L., Armanios E., Badir A. M. Theory of anisotropic thinwalled closedcross section beams. Composites Engineering, Vol. 2, Issue 57, 1992, p. 411432.

Choi S. C., Park. J. S., Kim J. H. Active damping of rotating composite thinwalled beam susing MFC actuators and PVDF sensors. Composite Structures, Vol. 76, 2006, p. 362374.

Zhong Y., Zhou L. Z. Basic theory of analytical mechanics. Beijing, Metallurgical Industry Press, 2006, (in Chinese).
About this article
The research is funded by the National Natural Science Foundation of China (Grant Nos. 10972124, 11272190) and Shandong Provincial Natural Science Foundation of China (Grant No. ZR2011EEM031).