Approximate Mode-based Simulation of Composite Wind Turbine Blade Vibrations using a Simplified Beam Model

N. Navadeh¹, I. O. Goroshko², Y. A. Zhuk² and A. S. Fallah^1,3,4,*

¹City and Guilds Building, Department of Aeronautics, South Kensington Campus, Imperial College London, London SW7 2AZ, UK

²Department of Theoretical and Applied Mechanics, Taras Shevchenko National University of Kyiv, Kiev 01601, Ukraine

³Howell Building, Department of Mechanical and Aerospace Engineering, Brunel University London, Uxbridge UB8 3PH, UK

⁴Institute of Computational Physics, Zürich University of Applied Sciences (ZHAW), 21 Wildbachstrasse, Winterthur 8401, Switzerland

E-mail: Arash.Soleiman-Fallah@brunel.ac.uk

* Corresponding Author

Received 04 February 2019; Accepted 16 August 2019; Publication 23 September 2019

Abstract

It is well-known that lower modes of vibration are responsible for a high percentage of the dynamic response. In this paper, the task of simulation of the dynamic response of the composite wind turbine blade on the basis numerical realisation of a developed low dimensional beam type model is considered. From the governing system of differential-algebraic equations of the simplified beam type model of the blade, and using the mode superposition approximation, the system of linear ordinary differential equations with respect to the coefficient functions of the modal representation was obtained. The developed program codes allow to simulate low frequency bending vibrations of wind turbine blades under different steady-state and transient loadings. The comparison of the simulation results obtained by the proposed simplified blade model with the results of the direct Finite Element Method (FEM) simulation shows their close agreement, which confirms the adequacy of the developed model and its mode-based approximation to the level of the requirements necessary in engineering practice. The presented approach to the creating low-dimensional simplified models of slender structures can therefore be useful in different fields of aerospace, civil, mechanical, and transport engineering.

Keywords: Low-dimensional beam model, transient response, vibration coupling, flapwise vibration, lead–lag vibration.

1 Introduction

Development of modern wind turbines and control of their dynamics in different operational modes need intensive computer simulation of stationary and transient vibrations of the turbine blades. Nowadays wind turbine blades are designed as long (tens of meters) slender composite shell structures with stiffeners, significantly inhomogeneous and non-prismatic (slenderising towards the free end (tip)) and with a variable angle of twist along the length. Modern FEM packages such as ANSYS, ABAQUS, or LS-DYNA allow to simulate mechanical deformation processes in such wind turbine blades with high accuracy (Chen and Chen, 2010), but in many cases, e.g. investigation of general structural dynamics, aero-elastic vibrations, preliminary design, and model based vibration control, such detailed costly models are redundant and can be replaced by simplified one-dimensional beam type models which allow to describe general blade deformations with appropriate accuracy and require much less simulation time or computational resources.

Parameters of the required simplified beam model are derived from results of physical or numerical (detailed FEM computations) experiments using different equalising and identification techniques. This approach is widely used for different aerospace structures (wings, fuselages, etc.) (Lee, 1995; Trivailo et al., 2006). But in the case of the wind turbine blades the situation is complicated due to considerable structural pre-twisting which does not render flap and lead-lag motions decomposable in the mathematical model.

In a previous article (Navadeh et al., 2017) an approach to construction of a beam type simplified model of a Horizontal Axis Wind Turbine (HAWT) composite blade, based on the results of more detailed FEM simulations of the blade, was proposed. In this method, the parameters of the model are obtained using the identification procedure from the FEM modal analysis data. This model allows effective description of low vibration bending modes of the blade taking into account the effects of coupling between flapwise and lead-lag vibrations. The present paper is devoted to simulation of transient vibrations of the wind turbine blade using the previously proposed simplified model, on the basis of the mode superposition method and evaluation of its effectiveness by comparison with the results of computations for the original shell type FEM model.

2 Low-Dimensional Simplified Beam-Mass Model of Turbine Blade

The real wind turbine blade is approximated by a piecewise homogeneous cantilever beam consisting of N weightless sections of length L_k with lumped masses m_k, $k=\overline{1,N}$, located at the connection points between sections and at the end point, as shown in Figure 1. To take into account the distribution of twist angle in the blade geometry along its length, principal axes of beam segments cross-sections x_k, y_k are rotated around longitudinal axis z by angles α_k (see Figure 2). The y-axis is oriented in the flap direction and the x-axis is in lead-lag bending direction. The x and y axes of the global coordinate system are oriented in the turbine rotor plane (lead-lag) and orthogonal to it (flap) directions, respectively.

The deformation of each beam section in local coordinates is described by standard Euler-Bernoulli equations:

$$\begin{array}{ccc}{({\mathit{EI}}_y)}_{{}_k}\frac{{\partial }^4{u}_k}{\partial z^4}=0& \text{and}& {({\mathit{EI}}_x)}_{{}_k}\frac{{\partial }^4{v}_k}{\partial z^4}=0,\end{array}$$

(1)

where (EI_y)_k and (EI_x)_k are the corresponding bending stiffnesses. Since the beam sections are connected to masses and are parts of the dynamic system, displacement components depend not only on z, the spatial coordinate, but also on time t: u_k = u_k(z, t), v_k = v_k(z, t). Thus, we are dealing with a discrete parameter model in which lumped masses are connected by beam elements with distributed stiffnesses described by Equations (1).

Figure 1 Schematic structure of the simplified beam-mass blade model.

Figure 2 Orientation of the principal axes of beam segments cross-section x_k, y_k with respect to axes of the global coordinate system x, y.

In the beam elements displacements u_k, v_k are presented by general solutions of homogeneous differential Equations (1) which are polynomials of the third order. Using local coordinate ξ_k (measured from the beginning of the segment in z direction), they can be written as:

$$\begin{array}{l}{u}_k=a_0{}^{(k)}+a_1{}^{(k)}{\xi }_k+a_2{}^{(k)}{\xi }_k^2+a_3{}^{(k)}{\xi }_k^3\hfill \\ v_k=b_0{}^{(k)}+b_1{}^{(k)}{\xi }_k+b_2{}^{(k)}{\xi }_k^2+b_3{}^{(k)}{\xi }_k^3\hfill \end{array}$$

(2)

As the generalised coordinates for the components of displacements, the coefficients $a_j^{(k)}$ and $b_j^{(k)}$ in (2) also are functions of time.

The components of basis unit vectors of the rotated k-th local coordinate system in terms of reference global coordinate system (subscript “0”) are presented as (see, for example, Reddy, 2006):

$$\begin{array}{cc}{\mathbf{e}}_x^{(k)}={(\cos {\alpha }_k,\sin {\alpha }_k)}_0,& {\mathbf{e}}_y^{(k)}={(\sin {\alpha }_k,\cos {\alpha }_k)}_0.\end{array}$$

(3)

Then, for the components of displacement vector for the k-th beam section in the global coordinate system one has

$$\begin{array}{l}{\mathbf{u}}_k={\left(u_k,v_k\right)}_k=u_k{\mathbf{e}}_k^{(k)}+v_k{\mathbf{e}}_y^{(k)}\hfill \\ =u_k{\left(\cos {\alpha }_k,\sin {\alpha }_k\right)}_0+v_k{\left(-\sin {\alpha }_k,\cos {\alpha }_k\right)}_0\hfill \\ ={\left(u_k\cos {\alpha }_k-v_k\sin {\alpha }_k,u_k\sin {\alpha }_k+v_k\cos {\alpha }_k\right)}_0,\hfill \end{array}$$

(4)

that corresponds completely to the relation given in (Reddy, 2006). Such transformations are universal and affect the components of any other vectors, such as moments and forces.

At the points of connection between sections continuity conditions are held i.e. the components, in the global coordinate system, of transversal displacements and their derivatives with respect to z (rotations) as well as components of bending moments are continuous. The components of shear forces have jumps due to the applied external forces and inertial forces from the lumped masses. At the free end (tip point), zero moment and shear force–inertial force conditions hold. At the clamped point of the beam (z = 0), the cantilever beam essential boundary conditions hold which give for the coefficients of the beam solutions (2) $a_0^{(1)}=a_1^{(1)}=b_0^{(1)}=b_1^{(1)}=0$. The conditions for displacements, rotations and moments have the form of linear algebraic equations, whereas the dynamic force conditions are ordinary differential equations. Together with expressions for the displacement components at the end point in global coordinate system u_END, v_END these equations form linear differential-algebraic (DAE) system with respect to components of the vector of governing functions of the problem:

$$c(t)=\{{\left\{a{j}^k\right\}}_{k=0}^N\cdot {\left\{b{j}^k\right\}}_{k=0}^N\cdot u_{\mathit{END}}\cdot v_{\mathit{END}}\}.$$

(5)

In the previous article (Navadeh et al., 2017) only free vibrations were considered, and the components of external forces applied to the lumped masses were omitted, thus the governing system of equations was homogeneous and in matrix notations had the form:

$$Ac=B\ddot{c}$$

(6)

where A and B are square and rectangular matrices of constant coefficients, respectively. The detailed representation of the governing system of equations in expanded form is given in (Navadeh, 2017) and will not be repeated here.

To obtain natural frequencies and eigenmodes i.e. shapes of natural vibration modes for the simplified beam-mass blade model the solutions of (6) are represented in exponential form $c(t)=\overline{c}{e}^{\lambda t}$, which leads to solving the generalised eigenvalue problem (Gruber, 2014):

$$A\overline{c}={\lambda }^{2}B\overline{c},$$

(7)

which gives the spectrum of eigenvalues, ${\lambda }_k^2$, and corresponding eigenvectors c_k. The natural cyclic frequencies of vibrations are calculated by the formula $f_k=\sqrt{{\lambda }_k^2}/(2\pi )$, and components of the c_k vectors are used to represent shapes of the normal vibration modes by Equations (2).

Note that in (Navadeh et al., 2017) the differential equations of the governing system, which are the component expressions in global coordinates of the second Newton’s law for the first N − 1 lumped masses of the considered beam-mass model, were not solved with respect to the higher-order (second) time derivatives of the $a_0^{k+1}$, $b_0^{k+1}$, $k=\overline{1,N-1}$, which are the components of displacements in local coordinate systems at the connection points between beam sections (see Equation (2)), i.e. the displacement components for the corresponding k-th masses. It’s not convenient for the further analysis, and these equations should be transformed into a simpler form.

From the conditions of the force balance at connection points between beam segments $(k=\overline{1,N-1})$ (Equation (7), (8) in (Navadeh et al., 2017)) and taking into account components $F_x^{(k)}(t)$, $F_y^{(k)}(t)$ of external forces applied at these points, we have following dynamic equations (hereinafter dots denote temporal differentiation):

$$\begin{array}{l}{m}_k({\ddot{a}}_0^{(k+1)}\cos {\alpha }_{k+1}-{\ddot{b}}_0^{(k+1)}\sin {\alpha }_{k+1})\hfill \\ =F_x^{(k)}(t)+6a_3^{(k)}{({\mathit{EI}}_y)}_k\cos {\alpha }_k-6a_3^{(k+1)}{({\mathit{EI}}_y)}_{k+1}\cos {\alpha }_{k+1}\hfill \\ -6b_3^{(k)}{({\mathit{EI}}_x)}_k\sin {\alpha }_k+6b_3^{(x+1)}{({\mathit{EI}}_k)}_{k+1}\sin {\alpha }_{k+1},\hfill \end{array}$$

(8)

$$\begin{array}{l}{m}_k({\ddot{a}}_0^{(k+1)}\sin {\alpha }_{k+1}+{\ddot{b}}_0^{(k+1)}\cos {\alpha }_{k+1})\hfill \\ =F_y^{(k)}(t)+6a_3^{(k)}{({\mathit{EI}}_y)}_k\sin {\alpha }_k-6a_3^{(k+1)}{({\mathit{EI}}_y)}_{k+1}\sin {\alpha }_{k+1}\hfill \\ +6b_3^{(k)}{({\mathit{EI}}_x)}_k\cos {\alpha }_k-6b_3^{(k+1)}{({\mathit{EI}}_k)}_{k+1}\cos {\alpha }_{k+1}.\hfill \end{array}$$

(9)

Multiplying Equation (8) by cos α_k+1 and Equation (9) by sin α_k+1 and adding them, we obtain

$$\begin{array}{l}{m}_k{\ddot{a}}_0^{(k+1)}=F_x^{(k)}(t)\cos {\alpha }_{k+1}+F_y^{(k)}(t)\sin {\alpha }_{k+1}\hfill \\ +6a_3^{(k)}{({\mathit{EI}}_y)}_k\cos ({\alpha }_k-{\alpha }_{k+1})\hfill \\ -6a_3^{(k+1)}{({\mathit{EI}}_y)}_{k+1}-6b_3^{(k)}{({\mathit{EI}}_x)}_k\sin ({\alpha }_k-{\alpha }_{k+1}).\hfill \end{array}$$

(10)

Then, multiplying (9) by cos α_k+1 and (8) by sin α_k+1 and subtracting the second from the first, we obtain

$$\begin{array}{l}{m}_k{\ddot{b}}_0^{(k+1)}=F_y^{(k)}(t)\cos {\alpha }_{k+1}-F_x^{(k)}(t)\sin {\alpha }_{k+1}\hfill \\ +6a_3^{(k)}{(E{I}_y)}_k\sin ({\alpha }_k-{\alpha }_{k+1})+6b_3^{(k)}{(E{I}_x)}_k\cos ({\alpha }_k-{\alpha }_{k+1})\hfill \\ -6b_3^{(k+1)}{(E{I}_x)}_{k+1}.\hfill \end{array}$$

(11)

The dynamic equations at the end point of the beam have the form (see Equation (9) in (Navadeh et al., 2017))

$$m_N{\ddot{u}}_{END}=F_x^{(N)}(t)+6a_3^{(N)}{({\mathit{EI}}_y)}_N\cos {\alpha }_N-6b_3^{(N)}{({\mathit{EI}}_x)}_N\sin {\alpha }_N,$$

(12)

$$m_N{\ddot{v}}_{END}=F_y^{(N)}(t)+6a_3^{(N)}{({\mathit{EI}}_y)}_N\sin {\alpha }_N-6b_3^{(N)}{({\mathit{EI}}_x)}_N\cos {\alpha }_N.$$

(13)

Since in further consideration the mode superposition method is applied, we will not present here the sufficiently large algebraic equations which express in governing DAE system the continuity conditions (Navadeh et al., 2017), because the components of eigenvectors, used for the solution representation, satisfy these equations identically.

3 Mode Superposition Method for Simulation of Transient Dynamics of the Wind Turbine Blade

The considered simplified beam-mass model, the parameters of which were obtained using the identification procedure, can effectively represent low frequency bending vibrations of the blade with structural pre-twisting (Navadeh et al., 2017). Because of this, to investigate transient deformations of the blade using this model, instead of the direct solution of the linear DAE system, which describes the dynamic behaviour of the wind turbine blade, that needs utilization of special numerical methods (and moreover can give additional errors due to the model inaccuracy in the higher frequency range), it’s more effective to use the mode superposition technique.

In the framework of the mode superposition method (Bathe, 2014), the solution for the vector of governing functions c(t) (5) of the linear DAE system is represented as a linear combination of eigenvectors c_m, obtained from the generalized eigenvalue problem (6), with time-dependent coefficient functions d_m(t). To investigate low frequency bending vibrations of the blade we will use the superposition of only the first M modes in the range, where the using beam type model has sufficient agreement with more detailed FEM one:

$$c(t)\approx {\displaystyle \sum _{m=1}^M{d}_m(t)c_m.}$$

(14)

Substituting corresponding components of the representation (14) into dynamic equations (10)–(13), we obtain the following system of 2N ordinary differential equations with respect to the functions d_m(t):

$$\begin{array}{l}{m}_k{\displaystyle \sum _{m=1}^M{\ddot{d}}_m(t)a_{0m}^{(k+1)}}=F_x^{(k)}(t)\cos {\alpha }_{k+1}+F_y^{(k)}(t)\sin {\alpha }_{k+1}\hfill \\ +{\displaystyle \sum _{m=1}^M{d}_m(t)A_m^{(k)},}\hfill \\ m_N{\displaystyle \sum _{m=1}^M{\ddot{d}}_m(t)u_{ENDm}}=F_x^{(N)}(t)+{\displaystyle \sum _{m=1}^M{d}_m(t)A_m^{(N)},}\hfill \end{array}$$

(15)

$$\begin{array}{l}{m}_k{\displaystyle \sum _{m=1}^M{\ddot{d}}_m(t)b_{0m}^{(k+1)}}=F_y^{(k)}(t)\cos {\alpha }_{k+1}-F_x^{(k)}(t)\sin {\alpha }_{k+1}\hfill \\ +{\displaystyle \sum _{m=1}^M{d}_m(t)B_m^{(k)},}\hfill \\ m_N{\displaystyle \sum _{m=1}^M{\ddot{d}}_m(t)v_{ENDm}}=F_y^{(N)}(t)+{\displaystyle \sum _{m=1}^M{d}_m(t)B_m^N.}\hfill \end{array}$$

Here $k=\overline{1,N-1}$, $a_{0m}^{(k+1)}$, $b_{0m}^{(k+1)}$, u_ENDm, v_ENDm are components of the m-th eigenvector c_m, and

$$\begin{array}{l}{A}_m^{(k)}=6a_{3m}^{(k)}{({\mathit{EI}}_y)}_k\cos ({\alpha }_k-{\alpha }_{k+1})-6a_{3m}^{(k+1)}{({\mathit{EI}}_y)}_{k+1}\hfill \\ -6b_{3m}^{(k)}{({\mathit{EI}}_x)}_k\sin ({\alpha }_k-{\alpha }_{k+1}),\hfill \\ A_m^{(N)}=6a_{3m}^{(N)}{({\mathit{EI}}_y)}_N\cos {\alpha }_N-6b_{3m}^{(N)}{({\mathit{EI}}_x)}_N\sin {\alpha }_N,\hfill \end{array}$$

(16)

$$\begin{array}{l}{B}_m^{(k)}=6a_{3m}^{(k)}{({\mathit{EI}}_y)}_k\sin ({\alpha }_k-{\alpha }_{k+1})+6b_{3m}^{(k)}{({\mathit{EI}}_x)}_k\cos ({\alpha }_k\hfill \\ -{\alpha }_{k+1})-6b_{3m}^{(k+1)}{({\mathit{EI}}_x)}_{k+1},\hfill \\ B_m^{(N)}=6a_{3m}^{(N)}{({\mathit{EI}}_y)}_N\sin {\alpha }_N+6b_{3m}^{(N)}{({\mathit{EI}}_x)}_N\cos {\alpha }_N.\hfill \end{array}$$

Then, using the Galerkin projection technique, we multiply equations in (15) by corresponding components of j-th $(j=\overline{1,M})$ eigenvectors used for the solution approximation (14) and sum them. As the result we obtain M ordinary differential equations with respect to M coefficient functions d_m(t):

$$\begin{array}{l}{\displaystyle \sum _{k=1}^{N-1}{m}_k}{\displaystyle \sum _{m=1}^M{\ddot{d}}_m(t)a_{0m}^{(k+1)}{a}_{0j}^{(k+1)}+m_N}{\displaystyle \sum _{m=1}^M{\ddot{d}}_m(t)u_{\mathit{ENDm}}{u}_{\mathit{ENDj}}}\hfill \\ +{\displaystyle \sum _{k=1}^{N-1}{m}_k}{\displaystyle \sum _{m=1}^M{\ddot{d}}_m(t)b_{0m}^{(k+1)}{b}_{0j}^{(k+1)}+m_N}{\displaystyle \sum _{m=1}^M{\ddot{d}}_m(t)v_{\mathit{ENDm}}{v}_{\mathit{ENDj}}}\hfill \\ ={\displaystyle \sum _{m=1}^M{\ddot{d}}_m}(t)[{\displaystyle \sum _{k=1}^{N-1}{m}_k{a}_{0m}^{(k+1)}{a}_{0j}^{(k+1)}+m_N}{u}_{\mathit{ENDm}}{u}_{\mathit{ENDj}}\hfill \\ +{\displaystyle \sum _{k=1}^{N-1}{m}_k{b}_{0m}^{(k+1)}{b}_{0j}^{(k+1)}+m_N{v}_{\mathit{ENDm}}{v}_{\mathit{ENDj}}}]\hfill \\ ={\displaystyle \sum _{m=1}^M{d}_m}(t)[{\displaystyle \sum _{k=1}^{N-1}{A}_m^{(k)}{a}_{0j}^{(k+1)}+A_m^{(N)}{u}_{\mathit{ENDj}}}\hfill \\ +{\displaystyle \sum _{k=1}^{N-1}{B}_m^{(k)}{b}_{0j}^{(k+1)}+B_m^{(N)}{v}_{\mathit{ENDj}}}]\hfill \\ +{\displaystyle \sum _{k=1}^{N-1}\left(F_x^{(k)}(t)\cos {\alpha }_{k+1}+F_y^{(k)}(t)\sin {\alpha }_{k+1}\right)a_{0j}^{(k+1)}+F_x^N(t)u_{\mathit{ENDj}}}\hfill \\ +{\displaystyle \sum _{k=1}^{N-1}\left(F_y^{(k)}(t)\cos {\alpha }_{k+1}-F_x^{(k)}(t)\sin {\alpha }_{k+1}\right)b_{0j}^{(k+1)}+F_y^N(t)v_{\mathit{ENDj}}.}\hfill \end{array}$$

(17)

Due to the orthogonality of displacements in eigenvectors with mass weights we have:

$$\begin{array}{l}{\displaystyle \sum _{m=1}^M{\ddot{d}}_m(t)}[{\displaystyle \sum _{k=1}^{N-1}{m}_k{a}_{0m}^{(k+1)}{a}_{0j}^{(k+1)}+m_N{u}_{\mathit{ENDm}}{u}_{\mathit{ENDj}}}\hfill \\ +{\displaystyle \sum _{k=1}^{N-1}{m}_k{b}_{0m}^{(k+1)}{b}_{0j}^{(k+1)}+m_N{v}_{\mathit{ENDm}}{v}_{\mathit{ENDj}}}]\hfill \\ ={\ddot{d}}_j(t)[{\displaystyle \sum _{k=1}^{N-1}{m}_k{(a_{0j}^{(k+1)})}^2+m_N{u}_{\mathit{ENDj}}^2}\hfill \\ +{\displaystyle \sum _{k=1}^{N-1}{m}_k{(b_{0j}^{(k+1)})}^2+m_N{v}_{\mathit{ENDj}}^2}],\hfill \end{array}$$

(18)

that simplifies solving the ODE system (17), since its mass matrix is diagonal.

Finally, we can represent the obtained ODE system (17) in the matrix-vector form, convenient for further numerical analyses (Bathe, 2014):

$$\ddot{\mathbf{d}}={\mathbf{P}}^{-1}\mathbf{Qd}+{\mathbf{P}}^{-1}\mathbf{f}.$$

(19)

Here d is a column vector of d_j(t); and the diagonal components of the mass matrix P are given as follows;

$$P_j={\displaystyle \sum _{k=1}^{N-1}{m}_k{(a_{0j}^{(k+1)})}^2+m_N{u}_{\mathit{ENDj}}^2}+{\displaystyle \sum _{k=1}^{N-1}{m}_k{(b_{0j}^{(k+1)})}^2+m_N{v}_{\mathit{ENDj}}^2,}$$

(20)

i.e. P = diag(P_j), and ${\mathbf{P}}^{-1}=\mathit{diag}(P_j^{-1})$, the components of the stiffness matrix Q have the form

$$Q_{\mathit{jm}}={\displaystyle \sum _{k=1}^{N-1}{A}_m^{(k)}{a}_{0j}^{(k+1)}+A_m^{(N)}{u}_{\mathit{ENDj}}}+{\displaystyle \sum _{k=1}^{N-1}{B}_m^{(k)}{b}_{0j}^{(k+1)}+B_m^{(N)}{v}_{\mathit{ENDj}}},$$

(21)

and the components of the load vector f are expressed as

$$\begin{array}{l}{f}_j={\displaystyle \sum _{k=1}^{N-1}(F_x^{(k)}(t)\cos {\alpha }_{k+1}+F_y^{(k)}(t)\sin {\alpha }_{k+1})a_{0j}^{(k+1)}+F_x^N(t)u_{\mathit{ENDj}}}\hfill \\ +{\displaystyle \sum _{k=1}^{N-1}(F_y^{(k)}(t)\cos {\alpha }_{k+1}-F_x^{(k)}(t)\sin {\alpha }_{k+1})b_{0j}^{(k+1)}+F_y^N(t)v_{\mathit{ENDj}}.}\hfill \end{array}$$

(22)

4 Results of Numerical Simulation

For simulation of the transient dynamics of the wind turbine blade using the developed multi-segment beam-mass model in the frame described above i.e. mode superposition approximation, and for the graphical representation of the results of computations, the Matlab codes were developed. Also for setting loads in the FEM model of the blade, developed in ABAQUS (2014), and extracting simulation results the corresponding Python scripts were written.

To evaluate the effectiveness of the simplified model of the blade in comparison with the detailed FEM model we consider the vibrations of the blade under a step (in time) load in y direction, distributed along the top reference line of the FEM model, as shown in Figure 3. The magnitude of the load density varies linearly from the value −1000 N at the root of the blade to −200 N at its end. This distribution was selected as the simplest approximation to the wind load simulation that is based upon the data for the blade chords given in (Chen and Chen, 2010).

In the present study, the simulations on the simplified beam model were carried out for N number of segments, here N = 6, using the model parameters obtained through the identification procedure in the previous article (Navadeh et al., 2017). Coordinates of the end nodes z_k (m), values of bending stiffnesses EI_x and EI_y, twisting angles α_k (degrees) for the k-th sections, and lumped masses m_k (kg), located at nodes, are given in Table 1:

In mode superposition representation (14) the first four modes were used.

Figure 3 Load distribution along the blade.

Table 1 Beam model segment data

k	1	2	3	4	5	6
zk	9.648	14.251	21.7315	29.212	36.6935	44.175
(EIx)k	7.2953·109	4.0715·108	1.1509·108	2.3845·107	8.6304·106	6.7160·106
(EIy)k	1.3532·1010	5.8728·109	1.9179·109	6.2110·108	2.7655·108	1.6842·108
αk	−12.856	−10.263	−6.335	−3.665	−2.06	−1.52
mk	6985.084	2003.617	1.1509·108	2.3845·107	8.6304·106	6.7160·106

The transient dynamic analysis for the shell type model in ABAQUS was fulfilled by the method of direct time integration. (Explicit temporal integration)

The results of simulation for displacements (m) in the y (flap) direction at the nodes, where lumped masses are located, are shown on the time segment from 0 to 3 seconds in Figure 4 (graphs a–f correspond to consecutive k-th nodes, see Figure 1). Solid lines represent the displacements, obtained by the simplified model, whereas dashed lines represent the results of the direct FEM simulation.

The comparison of the graphs in Figure 4 demonstrates a close agreement of the results, obtained by the simulation of the beam dynamic response under a step in time impulse loading, using simplified beam-type blade model, with more precise FEM ones. This indicates that the utilized four-mode approximation is adequate and allows to describe low frequency vibrations of the blade with sufficient accuracy for practical purposes.

Besides, in Figure 5 the time dependencies of the coefficient functions d_m(t) of the mode superposition representation (14), which define the contributions of different modes into the solution, are presented as follows.

From Figure 5 it is seen that at the considered loading the magnitude of d₁(t) is almost twenty times greater than that of d₂(t), and much more than d₃(t) and d₄(t). The result of this, taking into account the shapes in y direction of the first and second normal modes (which are represented in (Navadeh, 2017) in Figure 5a, b), is domination of the first mode’s contribution to displacements at the fifth and sixth nodes (see Figures 4e and f, respectively). But at the fourth node, where the shape function of the first mode became near four times less at the beam end, when the magnitudes of the second and third mode shape functions have absolute values some more than a half of their maxima (note that these modes have similar shapes in y direction, but different frequencies), as can be seen in Figure 4d, the contribution of the last modes is much more significant. Lastly, considering the graph of displacements at the third node (Figure 4c) we see that the presence in the used four-mode approximation of the fourth mode allow to represent motion in time more detailed than it could be done if only three modes were taken into account.

It should be noted that at the second node the displacements, obtained by the simplified beam model, are about one and a half times greater in magnitude than the corresponding values computed using more detailed FEM model, but the shapes of their time dependencies remain similar. This is explained by a decrease in the accuracy of the classical beam model (1) which is used in the blade approximation. However this inaccuracy in many cases can be neglected, because the magnitude of vibrations observing here is ∼1% of the maximum magnitude at the blade’s end. The same can be said about the low relational accuracy of the solution by the beam approximated model at the first node, where its magnitude relative to the magnitude of vibrations at the blade’s end is ∼0.2% vise <0.04% after the FEM simulation.

Figure 4 Displacements (in meters) in flap direction at the nodes of the simplified beam model; time is shown in seconds.

Figure 5 Coefficient functions of the mode superposition representation: (a) d₁(t) (blue line), d₂(t) (red); (b) d₂(t) (blue), d₃(t) (red), d₄(t) (yellow).

5 Conclusions

Numerical study of transient dynamic response to pulse loading forms an important integral part of dynamical investigation required for analysis and design of wind turbine composite blade structures. In the present work the question of simulation of the wind turbine blade dynamics, on the basis numerical realisation of the approximated low dimensional beam type model (which was developed in (Navadeh et al., 2017)), is considered.

At first, the equations of motion, that are included in the governing DAE system of the simplified mass-beam model of the blade, were transformed to ODE’s, solved with respect to the higher-order time derivatives of the node displacement components of the vector of governing functions. This allows to simplify the application of the mode superposition technique for approximate simulation of transient wind turbine blade vibrations. This method was then applied using a truncated set of generalised coordinates.

In the framework of the mode superposition method, the solution for the vector of governing functions of the linear DAE system was represented as a linear combination of eigenvectors, obtained from the generalized eigenvalue problem, with time-dependent coefficient functions. It should be noted that in practice for the investigation of low frequency bending vibrations of the beam type structures frequently the approximations containing only the several first vibrational modes is used, and here the number of modes is limited in the range, where using beam type model lead to sufficient accuracy and reasonable agreement with the real blade behaviour (described more detailed by the FEM model). As the result of substitution of the mode superposition representation into governing DAE system and utilization of the Galerkin projection technique, the linear ODE system with respect to coefficient functions was obtained. Due to orthogonality of displacements in eigenvectors with mass weights this system has a diagonal mass matrix and can be solved by standard numerical methods.

For the evaluation of the accuracy of the developed simplified beam type model of the wind turbine blade with mode superposition approximation in comparison with the precise FEM model the vibrations of the blade under a temporally step loading, the spatial profile of which varies linearly along the blade, was considered.

The results of simulation of the transient dynamics of the blade in the frame described above as the approximate model approach show a close agreement with the results, obtained by the direct FEM simulation, which indicates the consistency of this approach and its adequacy for describing low frequency vibrations of wind turbine blades with sufficient accuracy for practical purposes.

Modelling and simulation approaches developed in (Navadeh et al., 2017) and in this paper can be applied not only for investigation of the wind turbine blade dynamics, but also in broader fields of aerospace, civil, mechanical, defense, and transport engineering (vibrations, stability and model based control of wings, fuselages, space launchers, tall buildings, bridges, pipelines and other one-dimensional structures).

References

ABAQUS. (2014). v. 6.14 User Guide.

Bathe, K. J. (2014). Finite Element Procedures, 2nd edition. Watertown, MA: K. J. Bathe.

Chen, K.-N. and Chen, P.-Y. (2010). Structural optimization of 3 MW wind turbine blades using a two-step procedure. International Journal for Simulation and Multidisciplinary Design Optimization, 4, 159–165.

Gruber, M. H. J. (2014). Matrix algebra for linear models. San Francisco, CA: Wiley.

Lee, U. (1995). Equivalent dynamic beam–rod models of aircraft wing structures. Aeronautical Journal, 99(990), 450–457.

Navadeh, N., Goroshko, I. O., Zhuk, Y.A. and Fallah, A. S. (2017). An FEM-based AI approach to model parameter identification for low vibration modes of wind turbine composite rotor blades, European Journal of Computational Mechanics, 26(5–6), 541–556.

Reddy, J.N. (2006). Theory and Analysis of Elastic Plates and Shells, 2nd edition, CRC Press.

Trivailo, P. M., Dulikravich, G. S., Sgarioto, D., and Gilbert, T. (2006). Inverse problem of aircraft structural parameter estimation: Application of neural networks. Inverse Problems in Science and Engineering, 14(4), 351–363.

Biographies

N. Navadeh is a final year Ph.D. researcher in Aeronautics Department of Imperial College London conducting research and specializing on mechanics of composite wind turbines blades. He completed his B.Sc. in Civil Engineering jointly at Amirkabir University (Tehran Polytechnic) and University of Birmingham in 2008. He then joined City University to study towards his M.Sc. degree in Structural engineering where he completed his M.Sc. in 2010. He is a member of the Royal Aeronautical Society (RAeS) and of American Institute of Aeronautics and Astronautics (AIAA).

I. O. Goroshko received his B.Sc. and M.Sc. Degrees in Mechanics from Taras Shevchenko National University of Kyiv, Ukraine; and Ph.D. in Mechanics and Mathematics from Timoshenko Institute of Mechanics, National Academy of Sciences of Ukraine. Since 2013, he has been working as a Senior Researcher at the Department of Theoretical and Applied Mechanics in his alma mater. His fields of interests are continuum mechanics, dynamics of solids, and finite element analysis with his special focus being the investigation of interactive modal vibration in solids and structures.

Y. A. Zhuk received his B.Sc. and M.Sc. degrees in Continuum Mechanics from Taras Shevchenko National University of Kyiv, Ukraine; and Ph.D. as well as D.Sc. in Mechanics and Mathematics from Timoshenko Institute of Mechanics, National Academy of Sciences of Ukraine. He received his Professorship in Mechanics and Mathematics in 2012 from Mykolaiv State University, Ukraine, and joined the Taras Shevchenko National University of Kyiv the same year. Since 2013 he has been the Head of Department of Theoretical and Applied Mechanics in his alma mater. Prof. Y. A. Zhuk is the recipient of the numerous Fellowships (Fellowship of the President of Ukraine, Fellowship of Ukrainian National Academy of Sciences) and international awards (such as Grant of International Association of Academies of Sciences (CIS), Royal Society/NATO Fellowship, and EPSRC grants). As for his international collaborative work, he worked at Hong Kong University of Science and Technology, Imperial College London, and University of Aberdeen. Prof. Zhuk is a member of the National Committee on Theoretical and Applied Mechanics of Ukraine and a recipient of the Award for Excellence in Teaching from the National Academy of Sciences of Ukraine (2017). Prof. Zhuk has authored numerous scientific papers in Ukrainian and international journals and has co-authored the textbook in experimental mechanics awarded with the Prize of the National Academy of Sciences of Ukraine for the Best Textbook (2016). His fields of interests are the coupled thermos-electro-mechanics, nonlinear dynamics, and composite materials.

A. S. Fallah received his BSc in Civil Engineering from Iran University of Science and Technology in 1998 and his MSc degree in Structural Engineering from University of Tehran in 2002. In 2006, he was awarded a PhD in Structures by Imperial College London. He then worked as a Post-doctoral Research Fellow at Imperial College and King’s College London from 2007 until 2015. In 2015 he was appointed a Lecturer (Assistant Professor) in Aerostructures at Imperial College’s Aeronautics Department, a post he held until 2017. He is currently a Senior Lecturer (Associate Professor) in Dynamics and Vibration in Mechanical Engineering Department of Brunel University London. He is a professional member of the Royal Aeronautical Society and a member of the Institute for Materials and Manufacturing. He holds visiting positions at Imperial College London, Kiev Polytechnic Institute, and Zurich University of Applied Sciences.

European Journal of Computational Mechanics, Vol. 28_4, 307–324.
doi: 10.13052/ejcm1958-5829.2842
© 2019 River Publishers

Abstract

1 Introduction

2 Low-Dimensional Simplified Beam-Mass Model of Turbine Blade

3 Mode Superposition Method for Simulation of Transient Dynamics of the Wind Turbine Blade

4 Results of Numerical Simulation

5 Conclusions

References

Biographies