Dynamic Analysis of Composite Wind Turbine Blades as Beams: An Analytical and Numerical Study

Abstract

This study focuses on the dynamic modelling and analysis of the wind turbine blades made of multiple layers of fibre reinforced composites and core materials. For this purpose, a novel three-dimensional analytical straight beam model for blades is formulated. This model assumes that the beam is made of functionally graded material (FGM) and has a variable and asymmetrical cross section. In this model, the blades are assumed to be thin, slender and long with a relatively straight axis. They have two main parts, namely the core and the shell. The so-called core consists of a lightweight isotropic foam material, which also adds significant damping to the system. The core material is covered by the shell, which is modelled using homogenous and orthotropic material assumptions as the structure is reinforced with continuous fibres. Therefore, the blades are modelled under a straight beam with varying cross-section assumptions, in which the effective elastic properties are acquired by homogenizing the cross section. The beam formulation for modelling the system is performed both analytically and numerically with the finite element method. The results of both methods are in well agreement. The maximum deviation between the results is found below 4%.

1. Introduction

The demand for energy increases day by day as the human population and industry grow fast. This continuous growth has brought the problem of sourcing the energy from the least environment-harming methods, in other words, sustainably. As one of the most popular sustainable methods, electrical energy production using kinetic energy of the wind is encouraged all around the world. Even though wind energy is relatively more sustainable, more environment-friendly than many of the commonly used energy generation methods, it may not be very cost-effective. For instance, the initial investment costs of wind turbines are generally very high. The energy generation outcomes and reliability issues of the wind turbines also make their feasibility challenging.

Larger blades and higher towers are required to increase the capacity of wind turbines [1]. However, the growing size of wind turbines makes the design process of these structures harder as the length/thickness ratio of the blades and the tower makes them too slender so that even their own weights become significant [2]. Thus, lightweight designs are found to be more effective in energy generation. To achieve such light and strong structures, composite materials usually become the best fits [1,3,4,5]. Blades usually are made of a strong and stiff shell as the skin and a stiffener web or lightweight core material. This structure reduces weight while preserving strength [1,6].

Among all components of the turbine structure, the blade is the most crucial part that requires extra care due to its slender geometry and the loads on itself [7]. Therefore, the modelling of the wind turbine blades is a major part of the design process of the blades. Early modelling studies usually include numerical methods such as the finite element method (FEM) since it can be considered the most general and complicated cases with acceptable approximations [8,9]. The loading combinations acting on the blade and consequently the stress-states that occur in the beam are quite complicated and vary along the beam axis. Thus, FEM becomes a very useful standardised method that can handle such complexities. Alongside these complexities, there are also aerodynamic/aeroelastic aspects of the actual blade. They usually limit the structural design of the blades and must be adapted accordingly to meet these complex requirements. Aerodynamics usually limits the outer form of the blade, and hence the structural improvements should be done with the materials rather than the geometry. Moreover, the major structural concern is often the fatigue issues, in order to ensure a reliable operation of wind turbines [6,10]. Therefore, wind loads should not excite the vibration modes of the blades to avoid excessive vibration displacements [11]. Briefly, the structural dynamic and aeroelastic behaviour of wind turbines should be considered during their design process [12,13,14].

With the consideration of the geometry of the blades, the beam or shell elements can be used instead of three-dimensional tetrahedral or hexahedral elements for predicting mechanical behaviour, as this can reduce the computational cost without affecting accuracy significantly [15]. On the other hand, shells are usually more representative of the outer stiff, layered structure of the blades [16]. For these reasons, shells are usually more accurate in the modelling of blades than beams since the shells represent the deformations of the cross section with better rightness [17]. The disadvantage of using shell models is their having more degrees of freedom than the beam models, making it computationally more expensive. To reduce the computational cost of using shells as well as the error that symmetric isotropic and homogeneous beams cause, more advanced beam formulations that are specifically designated to model certain structures, such as propellers and wind turbine blades by taking into account axial loading, shear deformations etc., are developed [16,18,19,20].

Apart from FEM, there are other numerical models using various improved beam theories for the representation of rotor blades that consider the variation of the cross-sectional, material properties and the loads acting on the structure [21,22,23]. As mentioned earlier, the most common beam modelling approaches consider only the cases with symmetrical cross sections and homogeneous and isotropic material. However, the blades have varying cross sections, and the composite materials used in the wind turbine blades usually consist of fibre reinforced polymers. These composites are strongly orthotropic with changing orientation angles of the reinforcing fibres along the beam axis. Such structures are called tow steered composites [24,25,26]. Although the wind turbine blades usually undergo large deformations, the most common beam formulations, namely the Euler–Bernoulli and Timoshenko beam theories for straight and uniform beams, only include small displacements. Under these assumptions, performing analytical solutions even for frames are possible [27]. To build accurate models of the wind turbine blades, more general formulations that can take into account the characteristics of the blades must be considered, both geometrically and material-wise [28,29]. It is possible to identify the complexities of wind turbine blades in two main groups as geometrical and material. The large displacements, the mechanical systems may undergo, as well as the arbitrary and varying cross sections, that the beams may have, can be included within the geometric complexities. On the other hand, changing material behaviour with the position and direction (functionally graded material (FGM) and orthotropy/anisotropy) belong to the material category.

Since the wind turbine blades are made of fibre-reinforced composites, they can behave as orthotropic beams [30,31]. Higher-order beam models are capable of representing the laminated orthotropic beams and besides, as they allow for analytical solutions for certain conditions [32]. Another way to solve orthotropic beams analytically is the elasticity approach [33]. When the analytical solutions are not performable, numerical methods such as the finite difference method can also be used for solving these problems [34]. Moreover, material properties vary depending not only on the direction but also on the position in a continuum, which can be expressed with a function formulated in space. As mentioned earlier, such materials are called FGM. Inside the beams, material properties may vary throughout the thickness and/or along the axis [35,36]. Even though the material properties are not constant, analytical solutions can be performed and are available in the literature [37,38].

The geometric complexities can originate from the large displacements that the blades display, as well as from the cross sections being asymmetrical and varying through the axis. Large displacements and rotations lead to nonlinear terms in the equations of motion of the beams. Thus, only some limited cases are analytically solvable [39,40]. These nonlinearities are also considered for the analytical solutions of composite beams in the literature [41,42,43]. Furthermore, varying cross sections decisively affect the stiffness of the whole beams, and arbitrary cross sections lead to coupling of different vibration and buckling modes (flexural and torsional) [44,45]. For the composite wind turbine blades, the shape of the cross section is complex, and the material properties vary, depending on the position throughout the thickness of the beam and along the axis of the beam. This of course influences the shear behaviour of the beam [46,47]. Multiple methods exist to extract cross-sectional properties [48,49,50,51].

After covering the literature on modelling of composite beams and blades, it can be claimed that there is limited research concerning beams that are made of functionally graded material (FGM) with a varying and asymmetrical cross section. In this study, a novel analytical three-dimensional straight beam model with a varying cross section made of axially FGM is formulated. Then, its modal analysis is done analytically. This beam model considers asymmetrical cross sections where $I_{xy}$ is nonzero. The cross-sectional properties, as well as the elastic properties, namely Young’s and shear moduli, are acquired using a three-dimensional finite element model and an open-source computer code, i.e., BECAS (Beam Cross-section Analysis Software), that can run on GNU Octave, also an open-source software [52,53]. Finally, the novel analytical formulation of the three-dimensional straight beam is validated with the modal analysis results of the three-dimensional finite element model.

2. Dynamic Modelling of Composite Wind Turbine Blades

Throughout this study, the DTU-10MW-RWT blade, which is shown in Figure 1, is selected as a reference model [54]. This blade is 86.3 m in length and is made of uniaxial, biaxial, and triaxial laminates with a balsa core.

In this section, both analytical and finite element models of the composite wind turbine blades are presented. Rotational speed effects were not considered in the modelling. Therefore, gyroscopic and centrifugal stiffening effects on the dynamic behaviour of the wind turbine blade are not included within the scope of this study.

2.1. Analytical Model

An analytical formulation for a straight beam with a varying cross section and made of functionally graded material is presented for modelling wind turbine blades. This formulation assumes that the beam has a perfectly straight axis and an asymmetrical, rigid and homogeneous cross section. Considering the geometry of the blade, these assumptions are not expected to introduce a significant amount of error on the natural frequencies.

To keep the beam model linear, the beam is divided into multiple regions, where the variations of cross-sectional and material properties can be neglected and the average values of these varying quantities are used for each region. Then, the system can be fully defined and solved analytically using the compatibility conditions between the regions.

The coordinate system used in the formulation is presented in Figure 2. Here, the $XYZ$ coordinate system is attached to $G$, the centre of mass (CM); the $X^{\prime }{Y}^{\prime }{Z}^{\prime }$ coordinate system is attached to $S$, the shear centre (SC); $e_1$ and $e_2$ are distances between $\left(x_{CM},y_{CM}\right)$, i.e., the coordinates of the centre of mass, and $\left(x_{SC},y_{SC}\right)$, coordinates of the shear centre.

In the analytical formulation, the equilibrium and compatibility equations in Reference [27] are modified for the beam, which has a nonuniform and asymmetrical cross section. The equilibrium equations are

(1)$\frac{d{F}_x}{dz}+q_x=0$

(2)$\frac{d{F}_y}{dz}+q_y=0$

(3)$\frac{d{F}_z}{dz}+q_z=0$

(4)$\frac{d{M}_x}{dz}-F_y+m_x=0$

(5)$\frac{d{M}_y}{dz}+F_x+m_y=0$

(6)$\frac{d{M}_Z}{dz}+m_Z=0$

where $q_x$, $q_y$ and $q_z$ are distributed forces; $m_x$, $m_y$ and $m_z$ are distributed moments; $F_x$, $F_y$ and $F_z$ are internal forces; and $M_x$, $M_y$ and $M_z$ are internal moments.

The compatibility equations for the beam can then be written as

(7)$\frac{du}{dz}-{\Omega }_y-{\gamma }_x=0$

(8)$\frac{dv}{dz}+{\Omega }_x-{\gamma }_y=0$

(9)$\frac{dw}{dz}-{\gamma }_z=0$

(10)$\frac{d{\Omega }_x}{dz}-{\omega }_x=0$

(11)$\frac{d{\Omega }_y}{dz}-{\omega }_y=0$

(12)$\frac{d{\Omega }_z}{dz}-{\omega }_z=0$

where $u$, $v$, and $w$ are translational displacements; ${\Omega }_x$, ${\Omega }_y$, and ${\Omega }_z$ are angular displacements; ${\omega }_x$, ${\omega }_y$, and ${\omega }_z$ are relative angles for a unit length of the beam; and ${\gamma }_x$, ${\gamma }_y$, and ${\gamma }_z$ are shear angles.

Under linear, elastic, and isotropic material assumptions, the constitutive equation can be written in the matrix form as follow:

(13)$\left[\begin{array}{c}{F}_x\\ F_y\\ F_z\end{array}\right]=\left[\begin{array}{ccc}\raisebox{1ex}{$GA$}\!\left/ \!\raisebox{-1ex}{$k_x$}\right.& 0& 0\\ 0& \raisebox{1ex}{$GA$}\!\left/ \!\raisebox{-1ex}{$k_y$}\right.& 0\\ 0& 0& EA\end{array}\right]\left[\begin{array}{c}{\gamma }_x\\ {\gamma }_y\\ {\gamma }_z\end{array}\right]$

(14)$\left[\begin{array}{c}{M}_x\\ M_y\\ M_z\end{array}\right]=\left[\begin{array}{ccc}E{I}_x& -E{I}_{xy}& 0\\ -E{I}_{xy}& E{I}_y& 0\\ 0& 0& G{I}_b\end{array}\right]\left[\begin{array}{c}{\omega }_x\\ {\omega }_y\\ {\omega }_z\end{array}\right]$

Here, $A$ is the area; $I_x$ and $I_y$ are the moment of inertias around $x$ and $y$ axes at the centre of mass (CM); $I_{xy}$ is the product-moment of inertia at CM; $I_b$ is the torsional constant; $E$ is the effective Young modulus; $G$ is the effective shear modulus; and $k_x$ and $k_y$ are shear correction factors around $x$ and $y$ axes, respectively.

Finally, by combining Equations (7)–(12) and Equations (13) and (14), the following equations are obtained as:

(15)$\frac{du}{dz}={\Omega }_y+\frac{F_x}{\raisebox{1ex}{$GA$}\!\left/ \!\raisebox{-1ex}{$k_x$}\right.}$

(16)$\frac{dv}{dz}=-{\Omega }_x+\frac{F_y}{\raisebox{1ex}{$GA$}\!\left/ \!\raisebox{-1ex}{$k_y$}\right.}$

(17)$\frac{dw}{dz}=\frac{F_z}{EA}$

(18)$\frac{d{\Omega }_x}{dz}=\frac{M_x{I}_y}{E\left(I_x{I}_y-I_{xy}^2\right)}+\frac{M_y{I}_{xy}}{E\left(I_x{I}_y-I_{xy}^2\right)}$

(19)$\frac{d{\Omega }_y}{dz}=\frac{M_x{I}_{xy}}{E\left(I_x{I}_y-I_{xy}^2\right)}+\frac{M_y{I}_x}{E\left(I_x{I}_y-I_{xy}^2\right)}$

(20)$\frac{d{\Omega }_z}{dz}=\frac{M_z}{G{I}_b}.$

To obtain the eigenvalue problem, D’Alembert’s principle is employed, and the distributed forces and distributed moments in Equations (1)–(6) are replaced by the inertia terms:

(21)$q_x=-\mu \frac{{\partial }^{2}u}{\partial t^2}+\mu e_2\frac{{\partial }^2{\Omega }_z}{\partial t^2}$

(22)$q_y=-\mu \frac{{\partial }^{2}v}{\partial t^2}-\mu e_1\frac{{\partial }^2{\Omega }_z}{\partial t^2}$

(23)$q_z=-\mu \frac{{\partial }^{2}w}{\partial t^2}$

(24)$m_x=-\frac{\mu }{A}{I}_x\frac{{\partial }^2{\Omega }_x}{\partial t^2}+\frac{\mu }{A}{I}_{xy}\frac{{\partial }^2{\Omega }_y}{\partial t^2}$

(25)$m_y=-\frac{\mu }{A}{I}_y\frac{{\partial }^2{\Omega }_y}{\partial t^2}+\frac{\mu }{A}{I}_{xy}\frac{{\partial }^2{\Omega }_x}{\partial t^2}$

(26)$m_Z=-\frac{\mu }{A}{I}_p\frac{{\partial }^2{\Omega }_z}{\partial t^2}-\mu e_1\frac{{\partial }^{2}v}{\partial t^2}+\mu e_2\frac{{\partial }^{2}u}{\partial t^2}$

Here, $I_p$ is the polar moment of inertia at CM, and $\mu$ is the mass per length. Under the harmonic vibration assumption, Equations (21)–(26) can be rewritten as

(27)$q_x=\mu {\omega }^{2}u-\mu e_2{\omega }^2{\Omega }_z$

(28)$q_y=\mu {\omega }^{2}v+\mu e_1{\omega }^2{\Omega }_z$

(29)$q_z=\mu {\omega }^{2}w$

(30)$m_x=\frac{\mu }{A}{I}_x{\omega }^2{\Omega }_x-\frac{\mu }{A}{I}_{xy}{\omega }^2{\Omega }_y$

(31)$m_y=\frac{\mu }{A}{I}_y{\omega }^2{\Omega }_y-\frac{\mu }{A}{I}_{xy}{\omega }^2{\Omega }_x$

(32)$m_Z=\frac{\mu }{A}{I}_p{\omega }^2{\Omega }_z+\mu e_1{\omega }^{2}v-\mu e_2{\omega }^{2}u$

where $\omega$ represents the angular frequency. By substituting Equations (27)–(32) into Equations (1)–(6) and writing them together with Equations (15)–(20) in matrix form, a system of differential equations is then obtained as

(33)$\frac{d\mathbf{y}}{dz}=\mathbf{Ay}$

where $\mathbf{y}$ is written as

(34)${\mathbf{y}}^T={\left[\begin{array}{cccccccccccc}u& v& w& {\Omega }_x& {\Omega }_y& {\Omega }_z& F_x& F_y& F_z& M_x& M_y& M_z\end{array}\right]}^T$

and matrix $\mathbf{A}$ is given in Appendix A. The exact solution of the system of differential equations in the Equation (33) can be obtained only with constant cross-sectional and material properties. However, the blade has varying properties as mentioned earlier. Therefore, the blade is divided into several regions, assuming each region possesses constant cross-sectional and material properties. The varying properties of the beam are averaged using the average value theorem as given in Equation (35).

(35)$\overline{f}=\frac{1}{z_{i+1}-z_i}{\int }_{zi}^{z_{i+1}}f\left(z\right)dz$

Here, $f\left(z\right)$ is the function to be averaged; $z_{i+1}$ and $z_i$ are the start and endpoints of the region $i$, respectively; and $\overline{f}$ is the average value of the function $f\left(z\right)$.

The exact solution of Equation (33) is known as

(36)$\mathbf{y}=e^{\mathbf{A}z}{\mathbf{y}}_o.$

Here, ${\mathbf{y}}_o$ is the initial values vector. Following an algebraic set of equations are obtained as

(37)$\mathbf{F}\left(\omega \right){\mathbf{y}}_o=0$

by applying boundary and continuity conditions. Finally, the determinant of the coefficient matrix $\mathbf{F}\left(\omega \right)$ in the Equation (37) must be equal to zero for nontrivial solutions. The natural frequencies of the structure are obtained.

2.2. Finite Element Model

The finite element method is a useful tool for the dynamic modelling of complex structures. In addition, computer models created by using this method can then be imported into other programs such as Abaqus (Dassault Systémes, Vélizy-Villacoublay, France) for the finite element analyses and BECAS for the computation of geometry/material properties. Here, to validate the abovementioned formulation and create an input for BECAS, the finite element model of the blade provided by [52] was used.

The finite element model of the composite wind turbine blade was created with shell elements in Abaqus (Dassault Systémes, Vélizy-Villacoublay, France). To create the mesh structure, the eight-node shell elements with reduced integration (S8R) were used. The mesh structure of the composite blade is shown in Figure 3.

The model consists of eleven different regions and for each region, various composite layups are created by defining the region, material, thickness, and orientation of each ply. As boundary condition, it is assumed that the blade is fixed at its root.

3. Results

3.1. Structural Properties of the Blade

In composite structures, the geometry of the cross section is related to the number of layers and the thickness of each layer. Therefore, cross-sectional properties such as the area, the area moment of inertia, and the torsional constant should be computed accurately by considering these dimensional properties of the composite structure. To analyse the anisotropic and arbitrary-shaped beam cross sections, the open-source software BECAS was employed [52]. For 27 different $z$ coordinates, the required cross-sectional and material properties were obtained from BECAS, and these properties are expressed as continuous functions by applying curve-fitting process in GNU Octave. Functions used in the curve-fitting process are given in Appendix B. This can be clearly seen from Figure 4, Figure 5 and Figure 6; the fitted curves and BECAS data are in good agreement with each other.

Figure 4 clearly shows that there is a downward trend along the beam axis in the $A,I_x,I_y,I_p$, and $I_b$ parameters, while $I_{xy}$ increases until a certain length, approximately at 20 m. Surprisingly, $I_{xy}$ starts to decrease after that point. In Figure 5, it should be noted that $k_x$ reaches the maximum value at the blade tip; the same location is a minimum point for $k_y$. In addition, there are significant fluctuations in the effective Young and shear moduli along the beam axis, which can be seen from Figure 6. It was noted that these values start from relatively high values and reach their minimum between 20 and 40 m. It was also observed that there is a significant decrease in the effective Young modulus at the blade tip. The mass per length ($\mu$) has almost the same tendency with the area ($A$) along the beam axis.

3.2. Modal Analysis

Analytically, the first eight vibration modes of the composite blade were obtained after applying the analytical formulation using the generated curve-fitted functions as the cross-sectional and material properties for 500 regions along the beam axis. Similarly, the finite element modal analysis of the composite beam was also performed using the finite element software Abaqus (Dassault Systémes, Vélizy-Villacoublay, France). The analytical and finite element modal analysis results are presented in Table 1.

Based on analytical and FEM modal analyses, the vibration modes of the composite blade were identified as either flapwise bending, edgewise bending, coupled flapwise-edgewise bending, or torsional. The flapwise bending mode can be expressed as the bending mode about $x$ axis. Similarly, the edgewise bending mode is the bending mode about $y$ axis. The coupled flapwise-edgewise bending mode can be considered as a combination of the flapwise and edgewise bending modes, and it occurs in 3D space. The torsional mode can also be defined as the twist-dominant mode along $z$ axis. The first eight mode shapes of the reference blade were visualised using the finite element software Abaqus (Dassault Systémes, Vélizy-Villacoublay, France), as seen in Figure 7.

As can be seen from Table 1, the results of both analytical and finite element methods are in excellent agreement with the maximum error of approximately 3%. It was observed that the natural frequencies obtained with the analytical method are slightly higher than the FEM ones.

Figure 8 shows the effects of the number of regions that the beam is separated into for solution. Figure 8a shows that after 25 regions there is no significant change in the calculated first natural frequency of the beam. Figure 8b represents the CPU time elapsed to run the solution process. It should be noted that the computation was performed for the first 20 natural frequencies without using parallelisation and ran on a computer with an 8 GB RAM and Intel (Santa Clara, CA, United States) Core i7-6500U CPU @ 2.50 GHz (4 CPUs), ~2.6 GHz. It can be confidently claimed that after 100 regions, a significant rise in computational time was observed. Combining these two plots, it can be stated that 50 regions would be sufficient to achieve a satisfactory accuracy-computational time balance for this case. It is crucial to stress that the necessary number of regions will vary depending on the functions that represent the changing properties of the cross section and the material which the beam/blade is made of. The Abaqus (Dassault Systémes, Vélizy-Villacoublay, France) model requires 184 s to solve the first 20 modes of this blade structure on the same machine, and to keep a better comparison. The FE model was solved also without being in parallelisation.

4. Discussion

The error/deviation between the results of the two methods can be explained with the assumptions of the beam formulation, which includes the homogenisation of the material properties, averaging the cross-sectional properties throughout the regions along the beam axis and rigid motion of the cross section during deformation of the beam. More specifically, the reason that the analytical beam model has a slightly stiffer characteristic than the 3D FE model lies in the rigid cross-section assumption of the beam. It is claimed that during the deformation of the beam, the cross section of the beam preserves its initial shape, or in other words, it remains rigid and only undergoes rigid body motion, translation, and rotation. This leads to a stiffer behaviour than a real blade where the skin and the stiffening web make a thin-walled structure that can actually deform. However, the error introduced by this assumption does not cause the results to divert significantly from the numerical model, which is considered to be more realistic and computationally more expensive.

Apart from the natural frequencies, some coupled mode shapes are obtained. Here, the flexural (bending) and torsional modes are coupled. This occurs due to the formulation of the beam theory, which includes an asymmetrical cross section. Having an asymmetrical cross section causes the elastic, shear centres, centroid and centre of mass to be non-coincident. Particularly, the centre of mass is located further away from the rest due to both asymmetrical cross section and varying material properties throughout the cross section. This causes a twist during a bending motion or bending during torsional motion. Theoretically, all the bending and torsional modes are coupled. However, in all modes among the first eight (except for modes 6 and 8), one of the bending modes or the torsional mode is quite dominant. Therefore, these modes are named after the dominant mode character.

5. Conclusions

The main motivation of this study is to represent wind turbine blades with a more simplistic and analytical model. Hence, this study presented a novel beam formulation for three-dimensional free vibrations of beams made of axially FGM, with a varying and asymmetrical cross section. To validate the beam model, a 3D FE wind turbine model was used. The functions that define the cross-sectional and material properties of the blade were calculated using an open-source software, BECAS, which provided the necessary data from discrete points of the blade. The formulation evaluated the beam by dividing the beam into sections along the axis. Each region was handled with constant material and cross-sectional properties by calculating the average values for each quantity along the region. This allowed us to keep the model linear and easy to solve.

Comparing the results of both models, it can be stated that both are consistent and in good agreement. The error was spotted as approximately 2–3% for each mode where the maximum error was obtained as roughly 3%. Furthermore, the coupling of torsional and bending modes were observed in the mode shapes obtained from the FEM which occurred due to the asymmetrical cross-sectional geometry and nonhomogeneous material properties.

Author Contributions

Conceptualization, M.T., A.T. and E.T.; methodology, M.T. and E.T.; software, M.T. and Ö.E.G.; validation, M.T., Ö.E.G. and E.T.; formal analysis, M.T. and Ö.E.G.; investigation, M.T. and Ö.E.G.; resources, A.T.; data curation, M.T. and Ö.E.G.; writing—original draft preparation, M.T., Ö.E.G., A.T. and E.T.; writing—review and editing, M.T., Ö.E.G., A.T. and E.T.; visualization, Ö.E.G.; supervision, E.T.; project administration, E.T. and A.T. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A

(A1)

$\mathbf{A}=\left[\begin{array}{cccccccccccc}\hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill \frac{k_x}{GA}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \frac{k_y}{GA}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \frac{1}{EA}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \frac{I_y}{E\left(I_x{I}_y-I_{xy}^2\right)}\hfill & \hfill \frac{I_y}{E\left(I_x{I}_y-I_{xy}^2\right)}\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \frac{I_y}{E\left(I_x{I}_y-I_{xy}^2\right)}\hfill & \hfill \frac{I_x}{E\left(I_x{I}_y-I_{xy}^2\right)}\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \frac{1}{G{I}_b}\hfill \\ \hfill -\mu {\lambda }^2\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \mu e_2{\lambda }^2\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill -\mu {\lambda }^2\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -\mu e_1{\lambda }^2\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill -\mu {\lambda }^2\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -\frac{\mu }{A}{I}_x{\lambda }^2\hfill & \hfill \frac{\mu }{A}{I}_{xy}{\lambda }^2\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \frac{\mu }{A}{I}_{xy}{\lambda }^2\hfill & \hfill -\frac{\mu }{A}{I}_y{\lambda }^2\hfill & \hfill 0\hfill & \hfill -1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill \mu e_2{\lambda }^2\hfill & \hfill -\mu e_1{\lambda }^2\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -\frac{\mu }{A}{I}_p{\lambda }^2\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \end{array}\right]$

Appendix B

(A2)

$\left\{A,I_b,I_p,I_{xx},I_{yy}\right\}=a_1{e}^{-{\left(\frac{z-b_1}{c_1}\right)}^2}+a_2{e}^{-{\left(\frac{z-b_2}{c_2}\right)}^2}$

(A3)

$\left\{\mu ,G\right\}=a_1{e}^{-{\left(\frac{z-b_1}{c_1}\right)}^2}+a_2{e}^{-{\left(\frac{z-b_2}{c_2}\right)}^2}+a_3{e}^{-{\left(\frac{z-b_3}{c_3}\right)}^2}$

(A4)

$\begin{array}{l}\left\{I_{xy},x_{CM}\right\}=a_0+a_1\cos zw+b_1\sin zw+a_2\cos 2zw+b_2\sin 2zw+a_3\cos 3zw\\ +b_3\sin 3zw+a_4\cos 4zw+b_4\sin 4zw+a_5\cos 5zw\\ +b_5\sin 5zw+a_6\cos 6zw+b_6\sin 6zw+a_7\cos 7zw\\ +7\sin 7zw\end{array}$

(A5)

$\begin{array}{l}{k}_x=a_0+a_1\cos zw+b_1\sin zw+a_2\cos 2zw+b_2\sin 2zw+a_3\cos 3zw\\ +b_3\sin 3zw+a_4\cos 4zw+b_4\sin 4zw+a_5\cos 5zw\\ +b_5\sin 5zw+a_6\cos 6zw+b_6\sin 6zw\end{array}$

(A6)

$$\begin{gathered} k_y=a_0+a_1\cos zw+b_1\sin zw+a_2\cos 2zw+b_2\sin 2zw+a_3\cos 3zw \\ +b_3\sin 3zw+a_4\cos 4zw+b_4\sin 4zw \end{gathered}$$

(A7)

$y_{CM}=a_1\sin \left(b_{1}z+c_1\right)+a_2\sin \left(b_{2}z+c_2\right)+a_3\sin \left(b_{3}z+c_3\right)+a_4\sin \left(b_{4}z+c_4\right)$

(A8)

$x_{SC}=a_1\sin \left(b_{1}z+c_1\right)+a_2\sin \left(b_{2}z+c_2\right)+a_3\sin \left(b_{3}z+c_3\right)$

(A9)

$$\begin{gathered} \left\{y_{SC},E\right\}=a_1\sin \left(b_{1}z+c_1\right)+a_2\sin \left(b_{2}z+c_2\right)+a_3\sin \left(b_{3}z+c_3\right) \\ +a_4\sin \left(b_{4}z+c_4\right)+a_5\sin \left(b_{5}z+c_5\right) \end{gathered}$$

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Table

Figure 1. DTU-10MW-RWT blade [54].

Figure 2. The coordinate system used in the formulation.

Figure 3. Geometry and mesh structure of the composite blade.

Figure 4. The variations of area and area moment of inertia properties of the blade.

Figure 4. The variations of area and area moment of inertia properties of the blade.

Figure 5. The variations of shear correction factors and coordinates of centre of mass (CM) and shear centre (SC) of the blade.

Figure 6. The variations of material properties of the blade.

Figure 7. Finite element mode shapes of the DTU-10MW-RWT blade.

Figure 8. Dependence on the number of regions on the beam: (a) the first natural frequency (b) elapsed CPU time.

Table 1

Vibration modes of DTU-10MW-RWT blade.

Mode	AnalyticalNat. Frequency [Hz]	FEMNat. Frequency [Hz]	Mode Type	Error [%]
1	0.63	0.61	Flapwise bending	3.17
2	0.98	0.95	Edgewise bending	3.06
3	1.76	1.76	Flapwise bending	0
4	2.88	2.84	Edgewise bending	1.39
5	3.66	3.60	Flapwise bending	1.64
6	5.81	5.70	Coupled flapwise-edgewise bending	1.89
7	5.89	5.73	Torsional	2.72
8	6.28	6.13	Coupled flapwise-edgewise bending	2.39