Determination of natural frequencies and mode shapes of a wind turbine rotor blade using Timoshenko beam elements

When simulating a wind turbine, the lowest eigenmodes of the rotor blades are usually used to describe their elastic deformation in the frame of a multi-body system. In this paper, a finite element beam model for the rotor blades is proposed which is based on the transfer matrix method. Both static and kinetic field matrices for the 3-D Timoshenko beam element are derived by the numerical integration of the differential equations of motion using a Runge–Kutta fourth-order procedure. In the general case, the beam reference axis is at an arbitrary location in the cross section. The inertia term in the motion differential equation is expressed using appropriate shape functions for the Timoshenko beam. The kinetic field matrix is built by numerical integration applied on the approximated inertia term. The beam element stiffness and mass matrices are calculated by simple matrix operations from both field matrices. The system stiffness and mass matrices of the rotor blade model are assembled in the usual finite element manner in a global coordinate system accounting for the structural twist angle and possible pre-bending. The program developed for the above-mentioned calculations and the final solution of the eigenvalue problem is accomplished using MuPAD, a symbolic math toolbox in MATLAB^®. The natural frequencies calculated using generic rotor blade data are compared with the results proposed from the FAST and ADAMS software.

1 Introduction

Vibration of an elastic system refers to a limited reciprocating motion of a particle or an object of the system. Wind turbines operate in an unsteady environment which results in a vibrating response (see Manwell et al., 2009). They consist of long slender structures (rotor blades and tower), which result in resonant frequencies that should be taken into account during the initial design and construction. When the excitation frequency of the vibrating system is near any natural frequency, an undesirable resonant state occurs with large amplitudes, which may lead to damage or even the collapse of the wind turbine or its components. The vibration response, especially of the rotor blades, depends on the stiffness which is a function of the materials used, the design and the size (see Jureczko et al., 2005). Therefore, the estimation of natural frequencies in the early design phase plays an important role in avoiding resonance.

The eigenmodes associated with the lowest natural frequencies are employed as shape functions to describe the elastic deformation of the rotor blade beam model in the frame of the usual simulation of the wind turbine as a multi-body system. Generally, the first two bending eigenmodes in each direction (flapwise and edgewise) and optionally the two additional torsional eigenmodes are used. The determination of the lowest eigenmodes with sufficient numerical accuracy is the first step with respect to the modal superposition applied to the deformational motion of the rotor blades.

Due to the geometrical complexity of the blade cross section profiles and the extended use of composite materials, the exact calculation of natural frequencies in the design process takes a considerable amount of time and financial expense with regards to the 3-D modelling of the rotor blade using CAD software; hence a simplified finite element beam model is necessary. A twisted rotor blade is simplified into a cantilever beam with a non-uniform cross section. The structural twist angle is implemented by changing the orientation of the principal axis of the blade cross section along the length of the blade.

In the finite element formulation of beams two linear beam theories are established, the Euler–Bernoulli beam model and the Timoshenko beam model. Although the Euler–Bernoulli beam theory is widely used, the Timoshenko beam theory is considered to be better as it incorporates the effects of transverse shear and the rotational inertia on the dynamic response of the beam (see Kaya, 2006). In the classical approach of the finite element formulation for the free vibration analysis of beams, the stiffness and mass matrices are derived using interpolation functions derived from second- and fourth-order Hermite polynomials. The stiffness matrix is derived from the following equation (see Wu, 2013):

$$\begin{array}{}\text{(1)}& \mathbf{K}\left(e\right)=\int {\mathbf{B}}^{\mathrm{T}}{\mathbf{D}}_{\mathrm{m}}\mathbf{B}\mathrm{d}v,\end{array}$$

where K(e) is the element stiffness matrix, B is the strain matrix and D_m is the elasticity matrix for the beam. The element mass matrix of the beam (see Wu, 2013) is derived using the following equation:

$$\begin{array}{}\text{(2)}& \mathbf{M}\left(e\right)=\int \mathit{\rho }{\mathbf{a}}^{\mathrm{T}}\mathbf{a}\mathrm{d}v,\end{array}$$

where M(e) is the element mass matrix, ρ is the mass density, v is the volume and a is the matrix of interpolation functions.

Using the above-mentioned standard relations and appropriate shape functions for the Euler–Bernoulli beam and the Timoshenko beam, the stiffness matrix and consistent mass matrix for the finite beam element can be derived. However, an alternative to this procedure, based on the transfer matrix method for the beam theory (see Graf and Vassilev, 2006:69–88 and Stanoev 2007) is developed in the present article. The element stiffness matrix can be derived on the basis of numerical integration of the first-order ordinary differential equation system for the differential beam element. The associated mass matrix can be developed by numerical integration of the inertia term in the differential equation of motion. The numerical integration results in static and kinetic field matrices, from which the element stiffness and mass matrices can be easily assembled.

In the present article, the above-mentioned procedure is used to determine the element stiffness and element mass matrix for the Timoshenko beam. The interpolation functions used for deriving the mass matrix are based on Hermite polynomials according to Bazoune and Khulief (2003). The system stiffness and mass matrices for the rotor blade are assembled in a global coordinate basis in the usual finite element manner. The numerical solution of the associated eigenvalue problem for the system without damping is computed using computer algebra software (in the frame of MATLAB^®).

In Sects. 2 and 3, the proposed method is described in detail; in Sect. 4 the method is applied on a rotor blade structure with 288 degrees of freedom (DOFs). The results for the natural frequencies and the corresponding eigenmodes are compared with the results calculated using the FAST and ADAMS software.

2 Theoretical framework for the 3-D Timoshenko beam

2.1 Kinematic relationships

The general assumptions in the linear beam theory are as follows:

• a.

  The beam reference (longitudinal) axis is at any arbitrary location of the cross section.

• b.

  The only stresses that occur are normal stresses σ_x and shear stresses τ_xy, τ_xz.

• c.

  Cross section planes, which are initially normal to the longitudinal axis, will remain plane after deformation.

The geometrical representation of the deformation state of a beam cross section is shown in the Fig. 1. The deformations u_p, v_p and w_p at a cross-sectional point P are determined by three displacement functions u(x), v(x) and w(x) and three cross-sectional rotation angles φ_x(x), φ_y(x) and φ_z(x) – all of them are a function of the beam axis coordinate x. The differential equation system is derived in accordance with Stanoev (2013):

Figure 1 Deformation at the point ($x,y,z$) (see Andersen and Nielsen, 2008).

As seen in Fig. 1, the displacement vector u_p can be expressed at any cross section point P as

$$\begin{array}{}\text{(3)}& {\mathit{u}}_p=\left[\begin{array}{c}{u}_p\left(x,y,z\right)\\ v_p\left(x,y,z\right)\\ w_p\left(x,y,z\right)\end{array}\right]=\left[\begin{array}{c}u\left(x\right)-y{\mathit{\phi }}_z\left(x\right)+z{\mathit{\phi }}_y\left(x\right)\\ v\left(x\right)-z{\mathit{\phi }}_x\left(x\right)\\ w\left(x\right)+y{\mathit{\phi }}_x\left(x\right)\end{array}\right]\end{array}$$

The three components of the strains occurring in the beam can be expressed as the derivatives of the displacement functions u_p, v_p and w_p. The axial strain ε_xx and the two shear strain components γ_xz and γ_xy are given by the following:

$$\begin{array}{}\text{(4a)}& {\displaystyle }& {\displaystyle }{\mathit{\epsilon }}_{xx}={\displaystyle \frac{\partial u_p}{\partial x}}=u^{\prime }-y{\mathit{\phi }}_z^{\prime }+z{\mathit{\phi }}_y^{\prime }\text{(4b)}& {\displaystyle }& {\displaystyle }{\mathit{\gamma }}_{xz}={\displaystyle \frac{\partial u_p}{\partial z}}+{\displaystyle \frac{\partial w_p}{\partial x}}=\underset{{\mathit{\gamma }}_y}{\underbrace{{\mathit{\phi }}_y+w^{\prime }}}+y{\mathit{\phi }}_x^{\prime }\text{(4c)}& {\displaystyle }& {\displaystyle }{\mathit{\gamma }}_{xy}={\displaystyle \frac{\partial u_p}{\partial y}}+{\displaystyle \frac{\partial v_p}{\partial x}}=\underset{{\mathit{\gamma }}_y}{\underbrace{-{\mathit{\phi }}_z+v^{\prime }}}-z{\mathit{\phi }}_x^{\prime },\end{array}$$

where γ_z and γ_y are the constant shear strains that are not neglected in Timoshenko beam theory.

$$\begin{array}{}\text{(5a)}& {\displaystyle }& {\displaystyle }{\mathit{\gamma }}_y=-{\mathit{\phi }}_z+v^{\prime }\text{(5b)}& {\displaystyle }& {\displaystyle }{\mathit{\gamma }}_z={\mathit{\phi }}_y+w^{\prime }\end{array}$$

2.2 Principle of virtual work for internal forces

The virtual work components for internal forces corresponding to normal stresses and shear stresses are given by the following:

$$\begin{array}{ll}{\displaystyle }-\mathit{\delta }{W}_i& {\displaystyle }=\underset{x}{\int }\left\{\mathit{\delta }{u}^{\prime }N+\mathit{\delta }{\mathit{\phi }}_z^{\prime }{M}_z+\mathit{\delta }{\mathit{\phi }}_y^{\prime }{M}_y\right\}\mathrm{d}x+\underset{x}{\int }\left\{\mathit{\delta }{\mathit{\gamma }}_z{Q}_z\right\}\mathrm{d}x\\ \text{(6)}& {\displaystyle }& {\displaystyle }+\underset{x}{\int }\left\{\mathit{\delta }{\mathit{\gamma }}_y{Q}_y\right\}\mathrm{d}x+\underset{x}{\int }\left\{\mathit{\delta }{{\mathit{\phi }}^{\prime }}_x{M}_{\mathrm{TP}}\right\}\mathrm{d}x,\end{array}$$

where N is the axial force, M_z and M_y are bending internal moments, Q_y and Q_z are the corresponding shear forces and M_TP is the St. Venant torsional moment.

With the introduction of the constitutive relations of Hooke's material law for the stresses corresponding to the internal forces in Eq. (6) and expressing the strains using Eqs. (4a)–(4c) and (5a)–(5b), the virtual work relationship is reformulated as

$$\begin{array}{ll}{\displaystyle }-\mathit{\delta }{W}_i& {\displaystyle }=\underset{x}{\int }\mathit{\{}E\left[\underset{A}{\underbrace{\left(\underset{A}{\int }\mathrm{d}A\right)}}\cdot u^{\prime }-\underset{A_y}{\underbrace{\left(\underset{A}{\int }y\mathrm{d}A\right)}}\cdot {\mathit{\phi }}_z^{\prime }+\underset{A_z}{\underbrace{\left(\underset{A}{\int }z\mathrm{d}A\right)}}\cdot {\mathit{\phi }}_y^{\prime }\right]\mathit{\delta }{u}^{\prime }\\ {\displaystyle }& {\displaystyle }-E\left[\underset{A_y}{\underbrace{\left(\underset{A}{\int }y\mathrm{d}A\right)}}\cdot u^{\prime }-\underset{A_{yy}}{\underbrace{\left(\underset{A}{\int }{y}^{\mathrm{2}}\mathrm{d}A\right)}}\cdot {\mathit{\phi }}_z^{\prime }+\underset{A_{yz}}{\underbrace{\left(\underset{A}{\int }yz\mathrm{d}A\right)}}\cdot {\mathit{\phi }}_y^{\prime }\right]\mathit{\delta }{\mathit{\phi }}_z^{\prime }\\ {\displaystyle }& {\displaystyle }+E\left[\underset{A_z}{\underbrace{\left(\underset{A}{\int }z\mathrm{d}A\right)}}\cdot u^{\prime }-\underset{A_{yz}}{\underbrace{\left(\underset{A}{\int }yz\mathrm{d}A\right)}}\cdot {\mathit{\phi }}_z^{\prime }+\underset{A_{zz}}{\underbrace{\left(\underset{A}{\int }{z}^{\mathrm{2}}\mathrm{d}A\right)}}\cdot {\mathit{\phi }}_y^{\prime }\right]\mathit{\delta }{\mathit{\phi }}_y^{\prime }\\ {\displaystyle }& {\displaystyle }+G\left[\begin{array}{l}{A}_{sz}\left(w^{\prime }+{\mathit{\phi }}_y\right)\mathit{\delta }{w}^{\prime }+A_{sz}(w^{\prime }+{\mathit{\phi }}_y)\mathit{\delta }{\mathit{\phi }}_y\\ +A_{sy}\left(v^{\prime }-{\mathit{\phi }}_z\right)\mathit{\delta }{v}^{\prime }-A_{sy}(v^{\prime }-{\mathit{\phi }}_z)\mathit{\delta }{\mathit{\phi }}_z\end{array}\right]\\ \text{(7)}& {\displaystyle }& {\displaystyle }+G\underset{I_{\mathrm{T}}}{\underbrace{\left[\underset{A}{\int }(z^{\mathrm{2}}+y^{\mathrm{2}})\mathrm{d}A\right]}}{\mathit{\phi }}_x^{\prime }\mathit{\delta }{\mathit{\phi }}_x^{\prime }\mathit{\}}\mathrm{d}x\end{array}$$

Here, $A_{sz}={\mathit{\alpha }}_{sz}\cdot A$ and $A_{sy}={\mathit{\alpha }}_{sy}\cdot A$ are the shear areas in the z and y directions respectively, α_sz, α_sy are the corresponding shear correction coefficients, A is the cross-sectional area, A_y is the static moment with respect to the z axis, A_z is the static moment with respect to the y axis, A_yy is the moment of inertia with respect to the z direction, A_zz is the moment of inertia with respect to the y direction, A_yz is the deviation moment of inertia and I_T is the torsional moment of inertia.

After the coefficient comparison in Eqs. (6) and (7) the internal forces corresponding to the normal stresses can be expressed by introducing the cross sectional stiffness matrix EA:

$$\begin{array}{ll}{\displaystyle }\left[\begin{array}{c}N\\ -M_z\\ M_y\end{array}\right]& {\displaystyle }=\underset{\mathbf{EA}}{\underbrace{\left[\begin{array}{ccc}EA& E{A}_y& E{A}_z\\ E{A}_y& E{A}_{yy}& E{A}_{yz}\\ E{A}_z& E{A}_{yz}& E{A}_{zz}\end{array}\right]}}\cdot \left[\begin{array}{c}{u}^{\prime }\\ -{\mathit{\phi }}_z^{\prime }\\ {\mathit{\phi }}_y^{\prime }\end{array}\right]\\ \text{(8)}& {\displaystyle }& {\displaystyle }\Rightarrow \left[\begin{array}{c}{u}^{\prime }\\ -{\mathit{\phi }}_z^{\prime }\\ {\mathit{\phi }}_y^{\prime }\end{array}\right]=(\mathbf{EA}{)}^{-\mathrm{1}}\left[\begin{array}{c}N\\ -M_z\\ M_y\end{array}\right]\end{array}$$

The shear stress components in Eqs. (6) and (7) lead to the following relations:

$$\begin{array}{}\text{(9)}& {\displaystyle }{\displaystyle }{M}_x=M_{\mathrm{TP}}=G{I}_{\mathrm{T}}{\mathit{\phi }}_x^{\prime }\end{array}$$

$$\begin{array}{}\text{(10a)}& {\displaystyle }& {\displaystyle }{Q}_z=G{A}_{sz}\left(w^{\prime }+{\mathit{\phi }}_y\right)\text{(10b)}& {\displaystyle }& {\displaystyle }{Q}_y=G{A}_{sy}\left(v^{\prime }-{\mathit{\phi }}_z\right)\end{array}$$

The relation described in Eq. (10a) and (10b) implies that the chosen reference axis coincides with the shear centre due to neglected shear–torsion coupling terms in Eq. (7).

For the special case where the cross section coordinate system coincides with the principal axes, the deformation relationship in Eq. (8) reduces to

$$\begin{array}{ll}{\displaystyle }\left[\begin{array}{c}N\\ -M_z\\ M_y\end{array}\right]& {\displaystyle }=\left[\begin{array}{ccc}EA& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& E{A}_{YY}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& E{A}_{zz}\end{array}\right]\cdot \left[\begin{array}{c}{u}^{\prime }\\ -{\mathit{\phi }}_z^{\prime }\\ {\mathit{\phi }}_y^{\prime }\end{array}\right]\\ \text{(11)}& {\displaystyle }& {\displaystyle }\Rightarrow \left[\begin{array}{c}{u}^{\prime }\\ -{\mathit{\phi }}_z^{\prime }\\ {\mathit{\phi }}_y^{\prime }\end{array}\right]=(\mathbf{EA}{)}^{-\mathrm{1}}\left[\begin{array}{c}N\\ -M_z\\ M_y\end{array}\right]\end{array}$$

Figure 2 The finite beam element: internal forces and nodal DOFs.

2.3 Differential equation system

The virtual work relation in Eq. (7) is rewritten for the case where the beam coordinate system coincides with the principal axis of the cross section – see Eq. (11).

$$\begin{array}{ll}{\displaystyle }-\mathit{\delta }{W}_i& {\displaystyle }=\underset{x}{\int }\mathit{\{}\left(EA{u}^{\prime }\right)\mathit{\delta }{u}^{\prime }+\left(E{A}_{yy}{\mathit{\phi }}_z^{\prime }\right)\mathit{\delta }{\mathit{\phi }}_z^{\prime }+\left(E{A}_{zz}{\mathit{\phi }}_y^{\prime }\right)\mathit{\delta }{\mathit{\phi }}_y^{\prime }\\ {\displaystyle }& {\displaystyle }+\left(G{I}_{\mathrm{T}}{\mathit{\phi }}_x^{\prime }\right)\mathit{\delta }{\mathit{\phi }}_x^{\prime }+G{A}_{sz}\left(w^{\prime }+{\mathit{\phi }}_y\right)(\mathit{\delta }{w}^{\prime }+\mathit{\delta }{\mathit{\phi }}_y)\\ \text{(12)}& {\displaystyle }& {\displaystyle }+G{A}_{sy}\left(v^{\prime }-{\mathit{\phi }}_z\right)\left(\mathit{\delta }{v}^{\prime }-\mathit{\delta }{\mathit{\phi }}_z\right)\mathit{\}}\mathrm{d}x\end{array}$$

After the partial integration of Eq. (12) the cross section deformation relations and differential equilibrium conditions for the Timoshenko beam element are compiled in a first-order differential equation system (see Eqs. 13, 14). For the Timoshenko beam with an arbitrary beam reference axis at any point on the cross section (see Eq. 8), the system of differential equations can be expressed in the following form:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle \frac{\mathrm{d}}{\mathrm{d}x}}\left[\begin{array}{c}u\\ v\\ w\\ {\mathit{\phi }}_x\\ {\mathit{\phi }}_y\\ {\mathit{\phi }}_z\\ N\\ Q_y\\ Q_z\\ M_x\\ M_y\\ M_z\end{array}\right]=\\ {\displaystyle }& {\displaystyle }\left[\begin{array}{cccccccccccc}\mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& S_{\mathrm{11}}& \mathrm{0}& \mathrm{0}& \mathrm{0}& S_{\mathrm{13}}& -S_{\mathrm{12}}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}& {\displaystyle \frac{\mathrm{1}}{G{A}_{sy}}}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& -\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& {\displaystyle \frac{\mathrm{1}}{G{A}_{sz}}}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& {\displaystyle \frac{\mathrm{1}}{G{I}_{\mathrm{T}}}}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& S_{\mathrm{31}}& \mathrm{0}& \mathrm{0}& \mathrm{0}& S_{\mathrm{33}}& -S_{\mathrm{32}}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& -S_{\mathrm{21}}& \mathrm{0}& \mathrm{0}& \mathrm{0}& -S_{\mathrm{23}}& S_{\mathrm{22}}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& -\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\end{array}\right]\\ \text{(13)}& {\displaystyle }& {\displaystyle }\cdot \left[\begin{array}{c}u\\ v\\ w\\ {\mathit{\phi }}_x\\ {\mathit{\phi }}_y\\ {\mathit{\phi }}_z\\ N\\ Q_y\\ Q_z\\ M_x\\ M_y\\ M_z\end{array}\right]+\left[\begin{array}{c}\mathrm{0}\\ \mathrm{0}\\ \mathrm{0}\\ \mathrm{0}\\ \mathrm{0}\\ \mathrm{0}\\ -p_x\\ -p_y\\ -p_z\\ -m_T\\ -m_y\\ -m_z\end{array}\right]\end{array}$$

The differential equation system for the Timoshenko beam with the beam reference axis on principal axes can be represented in the following matrix form:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle \frac{\mathrm{d}}{\mathrm{d}x}}\left[\begin{array}{c}u\\ v\\ w\\ {\mathit{\phi }}_x\\ {\mathit{\phi }}_y\\ {\mathit{\phi }}_z\\ N\\ Q_y\\ Q_z\\ M_x\\ M_y\\ M_z\end{array}\right]=\\ {\displaystyle }& {\displaystyle }\left[\begin{array}{cccccccccccc}\mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& {\displaystyle \frac{\mathrm{1}}{EA}}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}& {\displaystyle \frac{\mathrm{1}}{G{A}_{sy}}}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& -\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& {\displaystyle \frac{\mathrm{1}}{G{A}_{sz}}}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& {\displaystyle \frac{\mathrm{1}}{G{I}_{\mathrm{T}}}}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& {\displaystyle \frac{\mathrm{1}}{E{A}_{zz}}}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& {\displaystyle \frac{\mathrm{1}}{E{A}_{zz}}}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& -\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\end{array}\right]\\ \text{(14)}& {\displaystyle }& {\displaystyle }\cdot \left[\begin{array}{c}u\\ v\\ w\\ {\mathit{\phi }}_x\\ {\mathit{\phi }}_y\\ {\mathit{\phi }}_z\\ N\\ Q_y\\ Q_z\\ M_x\\ M_y\\ M_z\end{array}\right]+\left[\begin{array}{c}\mathrm{0}\\ \mathrm{0}\\ \mathrm{0}\\ \mathrm{0}\\ \mathrm{0}\\ \mathrm{0}\\ -p_x\\ -p_y\\ -p_z\\ -m_T\\ -m_y\\ -m_z\end{array}\right]\end{array}$$

The entries S_ij in Eq. (13) are determined by the inversion of the cross-sectional stiffness matrix in Eq. (8):

$$\begin{array}{ll}{\displaystyle }\left[\begin{array}{c}{u}^{\prime }\\ -{\mathit{\phi }}_z^{\prime }\\ {\mathit{\phi }}_y^{\prime }\end{array}\right]& {\displaystyle }=\underset{\mathbf{S}}{\underbrace{{\left(\left[\begin{array}{ccc}EA& E{A}_y& E{A}_z\\ E{A}_y& E{A}_{yy}& E{A}_{yz}\\ E{A}_z& E{A}_{yz}& E{A}_{zz}\end{array}\right]\right)}^{-\mathrm{1}}}}\cdot \left[\begin{array}{c}N\\ -M_z\\ M_y\end{array}\right]\\ \text{(15)}& {\displaystyle }& {\displaystyle }=\underset{\mathbf{S}}{\underbrace{\left[\begin{array}{ccc}{S}_{\mathrm{11}}& S_{\mathrm{12}}& S_{\mathrm{13}}\\ S_{\mathrm{21}}& S_{\mathrm{22}}& S_{\mathrm{23}}\\ S_{\mathrm{31}}& S_{\mathrm{32}}& S_{\mathrm{33}}\end{array}\right]}}\cdot \left[\begin{array}{c}N\\ -M_z\\ M_y\end{array}\right]\end{array}$$

3 Alternative finite element formulation

In the classical finite element formulation, the beam stiffness matrix and the consistent beam mass matrix are derived by developing an approach for the displacement functions through shape (interpolation) functions, which consist of first- and third-order Hermite polynomials. In this section, an alternative finite element procedure is presented, based on the numerical Runge–Kutta fourth-order integration of the differential motion equations. The integration of the static part (the coefficient matrix in Eqs. 13 and 14 respectively) leads to the static field matrix, while the integration of the inertia terms in the equation of motion Eq. (16) results in a kinetic field matrix. In the last step of the integration, using both field matrices, the respective element stiffness and the element mass matrices can be calculated using simple matrix operations.

3.1 The differential equations of motion

The differential equations of motion for the differential beam element (see Fig. 2) can be written in a matrix form according to Stanoev (2007) and Müller (2012) as follows:

$$\begin{array}{ll}{\displaystyle }\left[\begin{array}{c}{\mathit{z}}_{\mathrm{1},x}\\ {\mathit{z}}_{\mathrm{2},x}\end{array}\right]& {\displaystyle }=\underset{\mathbf{A}}{\underbrace{\left[\begin{array}{cc}{\mathbf{A}}_{\mathrm{11}}& {\mathbf{A}}_{\mathrm{12}}\\ {\mathbf{A}}_{\mathrm{21}}& {\mathbf{A}}_{\mathrm{22}}\end{array}\right]}}\cdot \left[\begin{array}{c}{\mathit{z}}_{\mathrm{1}}\\ {\mathit{z}}_{\mathrm{2}}\end{array}\right]+\left[\begin{array}{c}{\mathit{b}}_{\mathrm{1}}\\ {\mathit{b}}_{\mathrm{2}}\end{array}\right]\\ \text{(16)}& {\displaystyle }& {\displaystyle }+\left[\begin{array}{c}\mathbf{0}\\ \mathbf{m}\end{array}\right]\cdot {\ddot{\mathit{z}}}_{\mathrm{1}}\end{array}$$

In the above-mentioned equation, matrix A is the coefficient matrix (see Eqs. 13, 14) and vector $\mathit{z}=\left[\begin{array}{c}{\mathit{z}}_{\mathrm{1}}\\ {\mathit{z}}_{\mathrm{2}}\end{array}\right]$ is the state vector, where

$$\begin{array}{ll}{\displaystyle }{\mathit{z}}_{\mathrm{1}}& {\displaystyle }={\left[\begin{array}{cccccc}u\left(x\right)& v\left(x\right)& w\left(x\right)& {\mathit{\phi }}_x\left(x\right)& {\mathit{\phi }}_y\left(x\right)& {\mathit{\phi }}_z\left(x\right)\end{array}\right]}^{\mathrm{T}}\\ \text{(16a)}& {\displaystyle }& {\displaystyle }\mathrm{\text{is the vector of the displacement functions, and}}{\displaystyle }{\mathit{z}}_{\mathrm{2}}& {\displaystyle }={\left[\begin{array}{cccccc}N\left(x\right)& Q_y\left(x\right)& Q_z\left(x\right)& M_x\left(x\right)& M_y\left(x\right)& M_z\left(x\right)\end{array}\right]}^{\mathrm{T}}\\ \text{(16b)}& {\displaystyle }& {\displaystyle }\mathrm{\text{is the vector of internal (section) force functions.}}\end{array}$$

The vector $\mathit{b}=\left[\begin{array}{c}{\mathit{b}}_{\mathrm{1}}\\ {\mathit{b}}_{\mathrm{2}}\end{array}\right]$ contains the known excitation forces; however, for an eigenvalue problem b=0. The coefficient matrix A and the state vector z in addition to excitation force b constitute the static part of the motion equation. The kinetic part of the motion equation can be expressed using a matrix of the interpolation functions and nodal acceleration vector as follows:

$$\begin{array}{}\text{(17)}& \left[\begin{array}{c}\mathbf{0}\\ \mathbf{m}\end{array}\right]\cdot {\ddot{\mathit{z}}}_{\mathrm{1}}=\underset{{\mathbf{b}}_{\mathrm{m}}}{\underbrace{\left[\begin{array}{c}\mathbf{0}\\ \mathbf{m}\end{array}\right]\left[{\mathbf{\Phi }}_{\mathrm{1}}\left(x\right){\mathbf{\Phi }}_{\mathrm{2}}\left(x\right)\right]}}\cdot \underset{{\ddot{\mathit{V}}}^R(e,t)}{\underbrace{\left[\begin{array}{c}\ddot{\mathit{V}}\left(a\right)\\ \ddot{\mathit{V}}\left(b\right)\end{array}\right]}},\end{array}$$

where ${\ddot{\mathit{z}}}_{\mathrm{1}}={\left[\begin{array}{cccccc}\ddot{u}& \ddot{v}& \ddot{w}& {\ddot{\mathit{\phi }}}_x& {\ddot{\mathit{\phi }}}_y& {\ddot{\mathit{\phi }}}_z\end{array}\right]}^{\mathrm{T}}$ is the vector of accelerations, [Φ₁(x)  Φ₂(x)]ϵR^(6×12) is the matrix of interpolation functions (see Sect. 3.3), $\ddot{\mathit{V}}\left(a\right),\ddot{\mathit{V}}\left(b\right)\mathit{ϵ}{R}^{(\mathrm{6}\times \mathrm{1})}$ is the vector with nodal accelerations and mϵR^(6×6) is the inertia matrix of the differential beam element (see Sect. 3.2).

3.2 The inertia matrix term

The inertia matrix in Eq. (18) implies that the distributed mass μ(x) (kg m⁻¹) is applied eccentrically at any location (y,z) in the cross section. The inertia matrix is expressed according to Stanoev (2013):

$$\begin{array}{}\text{(18)}& \mathbf{m}=\mathit{\mu }\cdot \left[\begin{array}{cccccc}\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& z& -y\\ \mathrm{0}& \mathrm{1}& \mathrm{0}& -z& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}& y& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& -z& y& (y^{\mathrm{2}}+z^{\mathrm{2}}+{\displaystyle \frac{{\mathrm{\Theta }}_p}{\mathit{\mu }}})& \mathrm{0}& \mathrm{0}\\ z& \mathrm{0}& \mathrm{0}& \mathrm{0}& (z^{\mathrm{2}}+{\displaystyle \frac{{\mathrm{\Theta }}_y}{\mathit{\mu }}})& -yz\\ -y& \mathrm{0}& \mathrm{0}& \mathrm{0}& -yz& (y^{\mathrm{2}}+{\displaystyle \frac{{\mathrm{\Theta }}_z}{\mathit{\mu }}})\end{array}\right],\end{array}$$

where Θ_p, Θ_y and Θ_z in (kg m) are the mass moments of inertia for the cross section.

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\mathrm{\Theta }}_y={\displaystyle \frac{\mathit{\mu }\cdot I_y}{A}}={\displaystyle \frac{\mathit{\mu }\cdot A_{zz}}{A}},{\mathrm{\Theta }}_z={\displaystyle \frac{\mathit{\mu }\cdot I_z}{A}}={\displaystyle \frac{\mathit{\mu }\cdot A_{yy}}{A}},\\ \text{(19)}& {\displaystyle }& {\displaystyle }{\mathrm{\Theta }}_p={\mathrm{\Theta }}_y+{\mathrm{\Theta }}_z\end{array}$$

3.3 Shape functions for the Timoshenko beam element

The acceleration term ${\ddot{\mathit{z}}}_{\mathrm{1}}$ in Eq. (17) is expressed using the product of Hermite interpolating polynomials and the nodal acceleration vectors $\ddot{\mathit{V}}\left(a\right)$ and $\ddot{\mathit{V}}\left(b\right)$.

Shape functions for axial and torsional deformations u(ξ) and φ_x(ξ), respectively, are derived using first-order polynomial as follows:

$$\begin{array}{}\text{(20)}& u\left(\mathit{\xi }\right)=a_{\mathrm{1}}+a_{\mathrm{2}}\mathit{\xi }=\underset{{\mathit{N}}_u}{\underbrace{\left[\mathrm{1}\mathit{\xi }\right]}}\underset{\mathit{a}}{\underbrace{\left[\begin{array}{c}{a}_{\mathrm{1}}\\ a_{\mathrm{2}}\end{array}\right]}}={\mathit{N}}_u\cdot \mathit{a},\mathit{\xi }={\displaystyle \frac{x}{L}}\end{array}$$

To express the coefficients a_j in terms of the nodal displacements $u\left(\mathit{\xi }=\mathrm{0}\right)=u_a$, $u\left(\mathit{\xi }=\mathrm{1}\right)=u_b$ or in terms of the torsion ${\mathit{\phi }}_x\left(\mathit{\xi }=\mathrm{0}\right)={\mathit{\phi }}_{x_a},{\mathit{\phi }}_x\left(\mathit{\xi }=\mathrm{1}\right)={\mathit{\phi }}_{x_b}$ the following relations are applied to Eq. (20):

$$\begin{array}{ll}{\displaystyle }\underset{\mathit{u}}{\underbrace{\left[\begin{array}{c}{u}_a\\ u_b\end{array}\right]}}& {\displaystyle }=\underset{\mathbf{S}}{\underbrace{\left[\begin{array}{cc}\mathrm{1}& \mathrm{0}\\ \mathrm{1}& \mathrm{1}\end{array}\right]}}\left[\begin{array}{c}{a}_{\mathrm{1}}\\ a_{\mathrm{2}}\end{array}\right]\\ \text{(21)}& {\displaystyle }& {\displaystyle }\to \left[\begin{array}{c}{a}_{\mathrm{1}}\\ a_{\mathrm{2}}\end{array}\right]={\mathbf{S}}^{-\mathrm{1}}\cdot \mathit{u}=\left[\begin{array}{cc}\mathrm{1}& \mathrm{0}\\ -\mathrm{1}& \mathrm{1}\end{array}\right]\left[\begin{array}{c}{u}_a\\ u_b\end{array}\right]\end{array}$$

Substituting Eq. (21) into Eq. (20) results in the shape function for axial deformation

$$\begin{array}{ll}{\displaystyle }u\left(\mathit{\xi }\right)& {\displaystyle }=\left[\mathrm{1}\mathit{\xi }\right]\underset{{\mathbf{G}}_u}{\underbrace{\left[\begin{array}{cc}\mathrm{1}& \mathrm{0}\\ -\mathrm{1}& \mathrm{1}\end{array}\right]}}\underset{{\mathit{v}}_u}{\underbrace{\left[\begin{array}{c}{u}_a\\ u_b\end{array}\right]}}={\mathit{N}}_u\cdot {\mathbf{G}}_u\cdot {\mathit{v}}_u\\ \text{(22)}& {\displaystyle }& {\displaystyle }=\underset{\mathrm{1}-\mathit{\xi }}{\underbrace{H_{u_{\mathrm{1}}}}}{u}_a+\underset{\mathit{\xi }}{\underbrace{H_{u_{\mathrm{2}}}}}{u}_b\end{array}$$

and in the following function for torsional deformation φ_x:

$$\begin{array}{}\text{(23)}& {\mathit{\phi }}_x\left(\mathit{\xi }\right)=H_{u_{\mathrm{1}}}{\mathit{\phi }}_{x_a}+H_{u_{\mathrm{2}}}{\mathit{\phi }}_{x_b}\end{array}$$

Figure 3 Definition of the dimensionless coordinate (ξ) of a beam element.

The relationships in Eq. (10a), Eq. (11) and the corresponding part of Eq. (14) are a starting point to derive approximation functions for bending deformation in the xz plane:

$$\begin{array}{}\text{(24)}& {\displaystyle }{Q}_z=G{A}_{sz}\left(w^{\prime }+{\mathit{\phi }}_y\right)=M_y^{\prime }=E{A}_{zz}{\mathit{\phi }}_y^{\prime \prime }\end{array}$$

Using the above-mentioned relation the expression for w^′ is given by

$$\begin{array}{ll}{\displaystyle }{w}^{\prime }& {\displaystyle }=-{\mathit{\phi }}_y+{\displaystyle \frac{E{A}_{zz}}{G{A}_{sz}}}{\mathit{\phi }}_y^{\prime \prime }=-{\mathit{\phi }}_y+\underset{{\mathit{\gamma }}_y}{\underbrace{{\mathit{\eta }}_y{\displaystyle \frac{L^{\mathrm{2}}}{\mathrm{12}}}{\mathit{\phi }}_y^{\prime \prime }}}\\ \text{(25)}& {\displaystyle }& {\displaystyle }\leftarrow {\mathit{\eta }}_y={\displaystyle \frac{\mathrm{12}E{A}_{zz}}{G{A}_{sz}{L}^{\mathrm{2}}}}\end{array}$$

The translational deformation function w(ξ) is approximated by a cubic polynomial function:

$$\begin{array}{ll}{\displaystyle }w\left(\mathit{\xi }\right)& {\displaystyle }=c_{\mathrm{0}}+c_{\mathrm{1}}\mathit{\xi }+c_{\mathrm{2}}{\mathit{\xi }}^{\mathrm{2}}+c_{\mathrm{3}}{\mathit{\xi }}^{\mathrm{3}}\\ \text{(26)}& {\displaystyle }& {\displaystyle }=\underset{{\mathit{N}}_w}{\underbrace{\left[\mathrm{1}\mathit{\xi }{\mathit{\xi }}^{\mathrm{2}}{\mathit{\xi }}^{\mathrm{3}}\right]}}\underset{\mathit{c}}{\underbrace{\left[\begin{array}{c}{c}_{\mathrm{0}}\\ c_{\mathrm{1}}\\ c_{\mathrm{2}}\\ c_{\mathrm{3}}\end{array}\right]}}={\mathit{N}}_w\cdot \mathit{c}\end{array}$$

Using the constant shear strain relation in Eqs. (5b) and (26) the polynomial expression for constant shear strain can be deduced as follows:

$$\begin{array}{}\text{(27)}& {\displaystyle }{\mathit{\gamma }}_z={\mathit{\eta }}_y{\displaystyle \frac{L^{\mathrm{2}}}{\mathrm{12}}}{\mathit{\phi }}_y^{\prime \prime }\text{, where}{\mathit{\phi }}_y^{\prime \prime }=-w^{\prime \prime \prime }=-{\displaystyle \frac{\mathrm{6}{c}_{\mathrm{3}}}{L^{\mathrm{3}}}}\end{array}$$

By including Eqs. (27) and (26) in Eq. (25) the polynomial expression for φ_y(ξ) results in

$$\begin{array}{ll}{\displaystyle }{\mathit{\phi }}_y\left(\mathit{\xi }\right)& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{L}}\left[\mathrm{0}-\mathrm{1}-\mathrm{2}\mathit{\xi }-{\displaystyle \frac{{\mathit{\eta }}_y}{\mathrm{2}}}-\mathrm{3}{\mathit{\xi }}^{\mathrm{2}}\right]\left[\begin{array}{c}{c}_{\mathrm{0}}\\ c_{\mathrm{1}}\\ c_{\mathrm{2}}\\ c_{\mathrm{3}}\end{array}\right]\\ \text{(28)}& {\displaystyle }& {\displaystyle }={\mathit{N}}_{{\mathit{\phi }}_y}\cdot \mathit{c}\end{array}$$

To determine the coefficients c_j the following boundary conditions are applied:

$$\begin{array}{ll}{\displaystyle }{\mathit{v}}_w& {\displaystyle }=\left[\begin{array}{c}{w}_a\\ {\mathit{\phi }}_{ya}\\ w_b\\ {\mathit{\phi }}_{yb}\end{array}\right]=\left[\begin{array}{c}w(\mathit{\xi }=\mathrm{0})\\ {\mathit{\phi }}_y(\mathit{\xi }=\mathrm{0})\\ w(\mathit{\xi }=\mathrm{1})\\ {\mathit{\phi }}_y(\mathit{\xi }=\mathrm{1})\end{array}\right]\\ \text{(29)}& {\displaystyle }& {\displaystyle }=\underset{{\mathbf{H}}_w}{\underbrace{\left[\begin{array}{cccc}\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& -{\displaystyle \frac{\mathrm{1}}{L}}& \mathrm{0}& -{\displaystyle \frac{{\mathit{\eta }}_y}{\mathrm{2}L}}\\ \mathrm{1}& \mathrm{1}& \mathrm{1}& \mathrm{1}\\ \mathrm{0}& -{\displaystyle \frac{\mathrm{1}}{L}}& -{\displaystyle \frac{\mathrm{2}}{L}}& -{\displaystyle \frac{\mathrm{3}+\frac{{\mathit{\eta }}_y}{\mathrm{2}}}{L}}\end{array}\right]}}\cdot \left[\begin{array}{c}{c}_{\mathrm{0}}\\ c_{\mathrm{1}}\\ c_{\mathrm{2}}\\ c_{\mathrm{3}}\end{array}\right]={\mathbf{H}}_w\cdot \mathit{c}\end{array}$$

The inversion of Eq. (29) yields

$$\begin{array}{ll}{\displaystyle }\mathit{c}& {\displaystyle }={\left({\mathbf{H}}_w\right)}^{-\mathrm{1}}\cdot {\mathit{v}}_w={\displaystyle \frac{\mathrm{1}}{\mathrm{1}+{\mathit{\eta }}_y}}\\ {\displaystyle }& {\displaystyle }\cdot \underset{{\mathbf{G}}_w}{\underbrace{\left[\begin{array}{cccc}{\mathit{\eta }}_y+\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ -{\mathit{\eta }}_y& {\displaystyle \frac{-L({\mathit{\eta }}_y+\mathrm{2})}{\mathrm{2}}}& {\mathit{\eta }}_y& {\displaystyle \frac{{\mathit{\eta }}_y}{\mathrm{2}}}L\\ -\mathrm{3}& {\displaystyle \frac{L({\mathit{\eta }}_y+\mathrm{4})}{\mathrm{2}}}& \mathrm{3}& {\displaystyle \frac{-L({\mathit{\eta }}_y-\mathrm{2})}{\mathrm{2}}}\\ \mathrm{2}& -L& -\mathrm{2}& -L\end{array}\right]}}\left[\begin{array}{c}{w}_a\\ {\mathit{\phi }}_{y_a}\\ w_b\\ {\mathit{\phi }}_{y_b}\end{array}\right]\\ \text{(30)}& {\displaystyle }& {\displaystyle }={\mathbf{G}}_w\cdot {\mathit{v}}_w\end{array}$$

The interpolation functions for $w\left(x,y,z\right)$ and ${\mathit{\phi }}_y\left(x,y,z\right)$, Eqs. (26) and (28), can be expressed by employing Eq. (30):

$$\begin{array}{ll}{\displaystyle }w\left(\mathit{\xi }\right)& {\displaystyle }=\underset{{\mathit{N}}_w}{\underbrace{\left[\mathrm{1}\mathit{\xi }{\mathit{\xi }}^{\mathrm{2}}{\mathit{\xi }}^{\mathrm{3}}\right]}}\\ {\displaystyle }& {\displaystyle }\underset{{\mathbf{G}}_w}{\underbrace{{\displaystyle \frac{\mathrm{1}}{{\mathit{\eta }}_y+\mathrm{1}}}\left[\begin{array}{cccc}{\mathit{\eta }}_y+\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ -{\mathit{\eta }}_y& {\displaystyle \frac{-L({\mathit{\eta }}_y+\mathrm{2})}{\mathrm{2}}}& {\mathit{\eta }}_y& {\displaystyle \frac{{\mathit{\eta }}_y}{\mathrm{2}}}L\\ -\mathrm{3}& {\displaystyle \frac{L({\mathit{\eta }}_y+\mathrm{4})}{\mathrm{2}}}& \mathrm{3}& {\displaystyle \frac{-L({\mathit{\eta }}_y-\mathrm{2})}{\mathrm{2}}}\\ \mathrm{2}& -L& -\mathrm{2}& -L\end{array}\right]}}\\ {\displaystyle }& {\displaystyle }\underset{{\mathit{v}}_w}{\underbrace{\left[\begin{array}{c}{w}_a\\ {\mathit{\phi }}_{y_a}\\ w_b\\ {\mathit{\phi }}_{y_b}\end{array}\right]}}\\ \text{(31)}& {\displaystyle }& {\displaystyle }=H_{w_{\mathrm{1}}}{w}_a+H_{w_{\mathrm{2}}}{\mathit{\phi }}_{y_a}+H_{w_{\mathrm{3}}}{w}_b+H_{w_{\mathrm{4}}}{\mathit{\phi }}_{y_b},\end{array}$$

where the products of both matrices N_w and G_w are introduced as functions $H_{w_j}(j=\mathrm{1},\mathrm{\dots },\mathrm{4})$.

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\mathit{\phi }}_y\left(\mathit{\xi }\right)=\underset{{\mathit{N}}_{{\mathit{\phi }}_y}}{\underbrace{{\displaystyle \frac{\mathrm{1}}{L}}\left[\mathrm{0}-\mathrm{1}-\mathrm{2}\mathit{\xi }-{\displaystyle \frac{{\mathit{\eta }}_y}{\mathrm{2}}}-\mathrm{3}{\mathit{\xi }}^{\mathrm{2}}\right]}}{\displaystyle \frac{\mathrm{1}}{{\mathit{\eta }}_y+\mathrm{1}}}\\ {\displaystyle }& {\displaystyle }\left[\begin{array}{cccc}{\mathit{\eta }}_y+\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ -{\mathit{\eta }}_y& {\displaystyle \frac{-L({\mathit{\eta }}_y+\mathrm{2})}{\mathrm{2}}}& {\mathit{\eta }}_y& {\displaystyle \frac{{\mathit{\eta }}_y}{\mathrm{2}}}L\\ -\mathrm{3}& {\displaystyle \frac{L({\mathit{\eta }}_y+\mathrm{4})}{\mathrm{2}}}& \mathrm{3}& {\displaystyle \frac{-L({\mathit{\eta }}_y-\mathrm{2})}{\mathrm{2}}}\\ \mathrm{2}& -L& -\mathrm{2}& -L\end{array}\right]\cdot \left[\begin{array}{c}{w}_a\\ {\mathit{\phi }}_{y_a}\\ w_b\\ {\mathit{\phi }}_{y_b}\end{array}\right]\\ \text{(32)}& {\displaystyle }& {\displaystyle }=H_{{\mathit{\phi }}_{y_{\mathrm{1}}}}{w}_a+H_{{\mathit{\phi }}_{y_{\mathrm{2}}}}{\mathit{\phi }}_{y_a}+H_{{\mathit{\phi }}_{y_{\mathrm{3}}}}{w}_b+H_{{\mathit{\phi }}_{y_{\mathrm{4}}}}{\mathit{\phi }}_{y_b}\end{array}$$

In Eq. (32), the functions $H_{{\mathit{\phi }}_{y_j}}$ ($j=\mathrm{1},\mathrm{\dots },\mathrm{4}$) are introduced in an analogous manner.

A similar method is used for determining the approximation functions v(ξ) and φ_z(ξ) for bending deformation in the xy plane.

Starting with Eqs. (10b), (11) and (14), the following is obtained:

$$\begin{array}{ll}\text{(33)}& {\displaystyle }{Q}_y& {\displaystyle }=G{A}_{sy}\left(v^{\prime }-{\mathit{\phi }}_z\right)=-M_z^{\prime }=-E{A}_{yy}{\mathit{\phi }}_z^{\prime \prime }{\displaystyle }{v}^{\prime }& {\displaystyle }={\mathit{\phi }}_z-{\displaystyle \frac{E{A}_{yy}}{G{A}_{sy}}}{\mathit{\phi }}_z^{\prime \prime }={\mathit{\phi }}_z+{\mathit{\eta }}_z{\displaystyle \frac{L^{\mathrm{2}}}{\mathrm{12}}}{\mathit{\phi }}_z^{\prime \prime }\\ \text{(34)}& {\displaystyle }& {\displaystyle }\leftarrow {\mathit{\eta }}_z={\displaystyle \frac{\mathrm{12}E{A}_{yy}}{G{A}_{sy}{L}^{\mathrm{2}}}}\end{array}$$

The approximations analogous to Eqs. (31) and (32) can be derived as follows:

$$\begin{array}{}\text{(35)}& {\displaystyle }& {\displaystyle }v\left(\mathit{\xi }\right)=H_{v_{\mathrm{1}}}{v}_a+H_{v_{\mathrm{2}}}{\mathit{\phi }}_{z_a}+H_{v_{\mathrm{3}}}{v}_b+H_{v_{\mathrm{4}}}{\mathit{\phi }}_{z_b}\text{(36)}& {\displaystyle }& {\displaystyle }{\mathit{\phi }}_z\left(\mathit{\xi }\right)=H_{{\mathit{\phi }}_{z_{\mathrm{1}}}}{v}_a+H_{{\mathit{\phi }}_{z_{\mathrm{2}}}}{\mathit{\phi }}_{z_a}+H_{{\mathit{\phi }}_{z_{\mathrm{3}}}}{v}_b+H_{{\mathit{\phi }}_{z_{\mathrm{4}}}}{\mathit{\phi }}_{z_b}\end{array}$$

The H_∗j functions developed in Eqs. (31), (32), (35) and (36) are “static” shape functions for the Timoshenko beam. Supposing dependence on time alone for the nodal displacement vectors V(a) and V(b), the matrix of the interpolation functions [Φ₁(x)  Φ₂(x)] in the inertia term Eq. (17) can be developed using Eqs. (31), (32), (35) and (36) (see Kusuma Chandrashekhara, 2018):

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\mathbf{\Phi }}_{\mathrm{1}}\left(x\right)=\left[\begin{array}{cccccc}{H}_{u_{\mathrm{1}}}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& H_{v_{\mathrm{1}}}& \mathrm{0}& \mathrm{0}& \mathrm{0}& H_{v_{\mathrm{2}}}\\ \mathrm{0}& \mathrm{0}& H_{w_{\mathrm{1}}}& \mathrm{0}& H_{w_{\mathrm{2}}}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& H_{u_{\mathrm{1}}}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& H_{{\mathit{\phi }}_{y_{\mathrm{1}}}}& \mathrm{0}& H_{{\mathit{\phi }}_{y_{\mathrm{2}}}}& \mathrm{0}\\ \mathrm{0}& H_{{\mathit{\phi }}_{z_{\mathrm{1}}}}& \mathrm{0}& \mathrm{0}& \mathrm{0}& H_{{\mathit{\phi }}_{z_{\mathrm{2}}}}\end{array}\right],\\ \text{(37)}& {\displaystyle }& {\displaystyle }{\mathbf{\Phi }}_{\mathrm{2}}\left(x\right)=\left[\begin{array}{cccccc}{H}_{u_{\mathrm{2}}}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& H_{v_s}& \mathrm{0}& \mathrm{0}& \mathrm{0}& H_{v_{\mathrm{4}}}\\ \mathrm{0}& \mathrm{0}& H_{w_s}& \mathrm{0}& H_{w_{\mathrm{4}}}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& H_{u_{\mathrm{2}}}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& H_{{\mathit{\phi }}_{y_s}}& \mathrm{0}& H_{{\mathit{\phi }}_{y_{\mathrm{4}}}}& \mathrm{0}\\ \mathrm{0}& H_{{\mathit{\phi }}_{z_s}}& \mathrm{0}& \mathrm{0}& \mathrm{0}& H_{{\mathit{\phi }}_{z_{\mathrm{4}}}}\end{array}\right],\end{array}$$

3.4 Numerical integration

The special form of the numerical Runge–Kutta fourth-order integration method applied here is described in detail in Müller and Wolf (1975) and Schenk (2012). The integration operator is applied to the equations of motion in Eq. (16), i.e.

$$\begin{array}{ll}{\displaystyle }\left[\begin{array}{c}{\mathit{z}}_{\mathrm{1},x}\\ {\mathit{z}}_{\mathrm{2},x}\end{array}\right]& {\displaystyle }=\mathbf{A}\cdot \left[\begin{array}{c}{\mathit{z}}_{\mathrm{1}}\\ {\mathit{z}}_{\mathrm{2}}\end{array}\right]+\left[\begin{array}{c}{\mathit{b}}_{\mathrm{1}}\\ {\mathit{b}}_{\mathrm{2}}\end{array}\right]\\ \text{(38)}& {\displaystyle }& {\displaystyle }+\underset{{\mathbf{b}}_m}{\underbrace{\left[\begin{array}{c}\mathbf{0}\\ \mathbf{m}\end{array}\right]\left[{\mathbf{\Phi }}_{\mathrm{1}}\left(x\right){\mathbf{\Phi }}_{\mathrm{2}}\left(x\right)\right]}}\cdot {\ddot{\mathit{V}}}^R\left(e,t\right)\end{array}$$

In order to gain sufficient numerical precision the beam axis needs to be divided into at least 20 integration intervals. The integration operator transfers the known state vector at beginning of the integration interval to the end of the interval. The integration procedure starts with the state vector at the first node a, i.e. at location (x=0):

$$\begin{array}{}\text{(39)}& \underset{\mathit{z}\left(a\right)}{\underbrace{\left[\begin{array}{c}{\mathit{z}}_{\mathrm{1}}(x=\mathrm{0})\\ {\mathit{z}}_{\mathrm{2}}(x=\mathrm{0})\end{array}\right]}}=\left[\begin{array}{ccc}\mathrm{1}& & \\ & \mathrm{\cdots }& \\ & & \mathrm{1}\end{array}\right]\mathit{z}\left(a\right)\end{array}$$

The integration operator is subsequently applied to the evaluated coefficient matrix A at each interval by excluding the initial state vector z(a). The result are static field matrices F(x,a), multiplicatively linked to z(a) see Eq. (40). Therefore, each F(x,a) matrix “transfers” the state vector at location (x=0) to the end x of the integration interval considered (transfer matrix method). In the frame of the integration procedure the actual field matrix F(x,a) serves column wise as an initial vector for the next interval, and the components of the state vector z(a) remain excluded. The beam “load” vector $\left[\begin{array}{c}{\mathit{b}}_{\mathrm{1}}\\ {\mathit{b}}_{\mathrm{2}}\end{array}\right]$, Eq. (38), evaluated in the actual interval, yields the $\left[\begin{array}{c}{\mathit{\beta }}_{\mathrm{1}}\\ {\mathit{\beta }}_{\mathrm{2}}\\ \mathrm{1}\end{array}\right]$ column in the F(x,a) matrix after integration – Eq. (40).

Figure 4 Eccentrically applied mass m_E at the x_AE point of the beam in the 3-D case.

Figure 5 First flapwise eigenmode.

Figure 6 First edgewise eigenmode.

The numerical integration of the inertia term b_m in Eq. (38) is carried out column wise analogously to the “load” vector, by excluding the nodal accelerations ${\ddot{\mathit{V}}}^R\left(e,t\right)$ – the results are kinetic field matrices H(x,a) at the end of each integration interval (at location x) (see Eq. 40):

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\begin{array}{cc}\underset{}{\underbrace{\left[\begin{array}{c}{\mathit{z}}_{\mathrm{1}}\left(x\right)\\ {\mathit{z}}_{\mathrm{2}}\left(x\right)\\ \mathbf{1}\end{array}\right]}}=& \underset{}{\underbrace{\left[\begin{array}{ccc}{\mathit{\alpha }}_{\mathrm{11}}& {\mathit{\alpha }}_{\mathrm{12}}& {\mathit{\beta }}_{\mathrm{1}}\\ {\mathit{\alpha }}_{\mathrm{21}}& {\mathit{\alpha }}_{\mathrm{22}}& {\mathit{\beta }}_{\mathrm{2}}\\ \mathbf{0}& \mathbf{0}& \mathbf{1}\end{array}\right]}}\cdot \left[\begin{array}{c}{\mathit{z}}_{\mathrm{1}}\left(a\right)\\ {\mathit{z}}_{\mathrm{2}}\left(a\right)\\ \mathbf{1}\end{array}\right]\\ \mathit{z}\left(x\right)=& \mathbf{F}(x,a)\cdot \mathit{z}\left(a\right)\end{array}\\ \text{(40)}& {\displaystyle }& {\displaystyle }\begin{array}{c}+\underset{}{\underbrace{\left[\begin{array}{ccc}{\mathbf{H}}_{\mathrm{11}}& {\mathbf{H}}_{\mathrm{12}}& \mathbf{0}\\ {\mathbf{H}}_{\mathrm{21}}& {\mathbf{H}}_{\mathrm{22}}& \mathbf{0}\\ \mathbf{0}& \mathbf{0}& \mathbf{1}\end{array}\right]}}\cdot \left[\begin{array}{c}\ddot{\mathit{V}}\left(a\right)\\ \ddot{\mathit{V}}\left(b\right)\\ \mathbf{0}\end{array}\right]\\ +\mathbf{H}(x,a)\cdot {\ddot{\mathit{V}}}^R\left(e\right)\end{array}\end{array}$$

This type of numerical integration allows (slightly) varying values of the coefficients of the A matrix, the b_m inertia term and the b vector along the beam axis – i.e. all stiffness, mass and external load values of the beam element may vary. After the last integration step at the second node b, at location (x=L), static L(e) and kinetic H(e) field matrices are obtained:

$$\begin{array}{}\text{(41a)}& {\displaystyle }{\displaystyle }\left[\begin{array}{c}{\mathit{z}}_{\mathrm{1}}\left(b\right)\\ {\mathit{z}}_{\mathrm{2}}\left(b\right)\\ \mathrm{1}\end{array}\right]=\mathbf{L}\left(e\right)\cdot \left[\begin{array}{c}{\mathit{z}}_{\mathrm{1}}\left(a\right)\\ {\mathit{z}}_{\mathrm{2}}\left(a\right)\\ \mathrm{1}\end{array}\right]+\mathbf{H}\left(e\right)\cdot \left[\begin{array}{c}\ddot{\mathit{V}}\left(a\right)\\ \ddot{\mathit{V}}\left(b\right)\\ \mathrm{0}\end{array}\right]\end{array}$$

where

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\mathbf{L}\left(e\right)=\left[\begin{array}{ccc}{\mathbf{L}}_{\mathrm{11}}& {\mathbf{L}}_{\mathrm{12}}& {\mathit{f}}_{\mathrm{1}}\\ {\mathbf{L}}_{\mathrm{21}}& {\mathbf{L}}_{\mathrm{22}}& {\mathit{f}}_{\mathrm{2}}\\ \mathbf{0}& \mathbf{0}& \mathrm{1}\end{array}\right],\\ \text{(41b)}& {\displaystyle }& {\displaystyle }\mathbf{H}\left(e\right)=\left[\begin{array}{ccc}{\mathbf{H}}_{\mathrm{11}}& {\mathbf{H}}_{\mathrm{12}}& \mathbf{0}\\ {\mathbf{H}}_{\mathrm{21}}& {\mathbf{H}}_{\mathrm{22}}& \mathbf{0}\\ \mathbf{0}& \mathbf{0}& \mathrm{1}\end{array}\right]\end{array}$$

According to Eqs. (13) and (14) the state variable z₁ represents six-component displacement vectors for V(a) and V(b), respectively, and the state variable z₂ represents six-component internal forces vectors S(a) and S(b), at locations (x=0) and (x=L) respectively.

Figure 7 Second flapwise eigenmode.

3.5 The element stiffness and mass matrices

By solving the matrix equation in Eq. (41a) for S(a) and S(b), and accounting for the definitions of internal forces in finite element beam formulation (see Fig. 2),

$$\begin{array}{}\text{(42)}& \mathit{F}\left(a\right)=-\mathit{S}\left(a\right),\mathit{F}\left(b\right)=\mathit{S}\left(b\right),\end{array}$$

one can derive the element stiffness matrix K(e), the element mass matrix M(e) and the element forces and moments F⁰ using simple matrix operations as shown in Eq. (43). Then the beam element relationships for the internal forces can be formulated as

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\left[\begin{array}{c}\mathit{F}\left(a\right)\\ \mathit{F}\left(b\right)\end{array}\right]=\left[\begin{array}{c}{\mathit{F}}^{\mathrm{0}}\left(a\right)\\ {\mathit{F}}^{\mathrm{0}}\left(b\right)\end{array}\right]+\underset{\mathbf{K}\left(e\right)}{\underbrace{\left[\begin{array}{cc}{\mathbf{K}}_{aa}& {\mathbf{K}}_{ab}\\ {\mathbf{K}}_{ba}& {\mathbf{K}}_{bb}\end{array}\right]}}\cdot \left[\begin{array}{c}\mathit{V}\left(a\right)\\ \mathit{V}\left(b\right)\end{array}\right]\\ {\displaystyle }& {\displaystyle }+\underset{\mathbf{M}\left(e\right)}{\underbrace{\left[\begin{array}{cc}{\mathbf{M}}_{aa}& {\mathbf{M}}_{ab}\\ {\mathbf{M}}_{ba}& {\mathbf{M}}_{bb}\end{array}\right]}}\cdot \left[\begin{array}{c}\ddot{\mathit{V}}\left(a\right)\\ \ddot{\mathit{V}}\left(b\right)\end{array}\right]=\underset{{\mathbf{F}}^{\mathrm{0}}}{\underbrace{\left[\begin{array}{c}{\mathbf{L}}_{\mathrm{12}}^{-\mathrm{1}}\cdot {\mathit{f}}_{\mathrm{1}}\\ {\mathit{f}}_{\mathrm{2}}-{\mathbf{L}}_{\mathrm{22}}\cdot {\mathbf{L}}_{\mathrm{12}}^{-\mathrm{1}}\cdot {\mathit{f}}_{\mathrm{1}}\end{array}\right]}}\\ {\displaystyle }& {\displaystyle }+\underset{\mathbf{K}\left(e\right)}{\underbrace{\left[\begin{array}{cc}{\mathbf{L}}_{\mathrm{12}}^{-\mathrm{1}}\cdot {\mathbf{L}}_{\mathrm{11}}& -{\mathbf{L}}_{\mathrm{12}}^{-\mathrm{1}}\\ {\mathbf{L}}_{\mathrm{21}}-{\mathbf{L}}_{\mathrm{22}}\cdot {\mathbf{L}}_{\mathrm{12}}^{-\mathrm{1}}\cdot {\mathbf{L}}_{\mathrm{11}}& {\mathbf{L}}_{\mathrm{22}}\cdot {\mathbf{L}}_{\mathrm{12}}^{-\mathrm{1}}\end{array}\right]}}\cdot \left[\begin{array}{c}\mathit{V}\left(a\right)\\ \mathit{V}\left(b\right)\end{array}\right]\\ {\displaystyle }& {\displaystyle }+\underset{\mathbf{M}\left(e\right)}{\underbrace{\left[\begin{array}{cc}{\mathbf{L}}_{\mathrm{12}}^{-\mathrm{1}}\cdot {\mathbf{H}}_{\mathrm{11}}& -{\mathbf{L}}_{\mathrm{12}}^{-\mathrm{1}}\cdot {\mathbf{H}}_{\mathrm{12}}\\ {\mathbf{H}}_{\mathrm{21}}-{\mathbf{L}}_{\mathrm{22}}\cdot {\mathbf{L}}_{\mathrm{12}}^{-\mathrm{1}}\cdot {\mathbf{H}}_{\mathrm{11}}& {\mathbf{H}}_{\mathrm{22}}-{\mathbf{L}}_{\mathrm{22}}\cdot {\mathbf{L}}_{\mathrm{12}}^{-\mathrm{1}}\cdot {\mathbf{H}}_{\mathrm{12}}\end{array}\right]}}\\ \text{(43)}& {\displaystyle }& {\displaystyle }\cdot \left[\begin{array}{c}\ddot{\mathit{V}}\left(a\right)\\ \ddot{\mathit{V}}\left(b\right)\end{array}\right]\end{array}$$

Figure 8 Mixed flap/edgewise eigenmode.

3.6 Single masses at eccentric positions

The numerical integration according to Runge–Kutta, described in Sec. 3.4, offers the possibility of including single load or mass quantities within a beam element. Single eccentric masses can be taken into account at the integration interval boundaries. In the local coordinate system of the beam element, at a general position vector ${\mathit{x}}_{AE}={\left[y_E{z}_E\right]}^{\mathrm{T}}$, the eccentric mass and the vector representation of dynamic equilibrium, Eq. (44a)–(44b), is as shown in Fig. 4. The beam reference axis is at point A, and vector ${\ddot{\mathit{V}}}_E$ represents the acceleration vector at the point of application (x,x_AE). With the help of the dynamic equilibrium conditions Eq. (44a)–(44b), additional inertia forces and moments due to eccentric mass can be determined (see Li, 2016):

$$\begin{array}{}\text{(44a)}& {\displaystyle }& {\displaystyle }-{\mathit{N}}_L+{\mathit{N}}_R=m_E{\ddot{\mathit{V}}}_E\text{(44b)}& {\displaystyle }& {\displaystyle }{\mathit{M}}_L+{\mathit{M}}_R=\sum _{i=\mathrm{1}}^{\mathrm{3}}\left({\mathit{\theta }}_{Ei}{\ddot{\mathit{\phi }}}_i{\mathit{e}}_i\right)+\left({\mathit{x}}_{AE}\times m_E{\ddot{\mathit{V}}}_E\right),\end{array}$$

where N_L and N_R are internal force vectors on the left and right in differential proximity to the point at location x, M_L and M_R are internal moment vectors on the left and right in differential proximity to the point at location x (see Fig. 4), m_E is the eccentric single mass, Θ_Ei are the mass moments of inertia of the single mass and ${\ddot{\mathit{\phi }}}_i$ are the angular accelerations at location x.

Figure 9 First torsional eigenmode.

The additional inertia matrix ${\mathbf{M}}_E\left(m_E,{\mathrm{\Theta }}_{Ei},x,{\mathit{x}}_{AE}\right)$, derived from Eq. (44a)–(44b), is analogous to the inertia matrix due to distributed mass in Eq. (18). During the numerical integration within the beam element (see Sect. 3.4) an additional eccentric inertia term has to be added to the kinetic field matrix at the end x (the point of application of eccentric mass) of the corresponding integration interval – see Eqs. (40) and (45).

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\left[\begin{array}{c}{\mathit{z}}_{\mathrm{1}}\left(x\right)\\ {\mathit{z}}_{\mathrm{2}}\left(x\right)\\ \mathbf{1}\end{array}\right]=\left[\begin{array}{ccc}{\mathit{\alpha }}_{\mathrm{11}}& {\mathit{\alpha }}_{\mathrm{12}}& {\mathit{\beta }}_{\mathrm{1}}\\ {\mathit{\alpha }}_{\mathrm{21}}& {\mathit{\alpha }}_{\mathrm{22}}& {\mathit{\beta }}_{\mathrm{2}}\\ \mathbf{0}& \mathbf{0}& \mathbf{1}\end{array}\right]\cdot \left[\begin{array}{c}\begin{array}{c}{\mathit{z}}_{\mathrm{1}}\left(a\right)\\ {\mathit{z}}_{\mathrm{2}}\left(a\right)\\ \mathbf{1}\end{array}\end{array}\right]\\ {\displaystyle }& {\displaystyle }+\left[\left[\begin{array}{ccc}{\mathbf{H}}_{\mathrm{11}}& {\mathbf{H}}_{\mathrm{12}}& \mathbf{0}\\ {\mathbf{H}}_{\mathrm{21}}& {\mathbf{H}}_{\mathrm{22}}& \mathbf{0}\\ \mathbf{0}& \mathbf{0}& \mathbf{1}\end{array}\right]+\left[\begin{array}{ccc}\mathbf{0}& \mathbf{0}& \mathbf{0}\\ {\mathbf{M}}_E{\mathbf{\Phi }}_{\mathrm{1}}& {\mathbf{M}}_E{\mathbf{\Phi }}_{\mathrm{2}}& \mathbf{0}\\ \mathbf{0}& \mathbf{0}& \mathbf{1}\end{array}\right]\right]\\ \text{(45)}& {\displaystyle }& {\displaystyle }\cdot \left[\begin{array}{c}\ddot{\mathit{V}}\left(a\right)\\ \ddot{\mathit{V}}\left(b\right)\\ \mathbf{0}\end{array}\right]\end{array}$$

Single masses do not usually appear in a rotor blade model, but the same finite element may be used for the modelling of wind turbine towers. In this case single masses within a finite beam element could represent bolted ring flange connections or the mass of any equipment, such as lifts etc.

Table 1 First three calculated (bolded values) flapwise and edgewise eigenfrequencies.

Table 2 Comparison between Timoshenko and Bernoulli beams with three variants for shear stiffness values.

Figure 10 Third flapwise eigenmode.

Figure 11 Second torsional eigenmode.

4 The eigenvalue problem

The system matrices for a rotor blade beam model are assembled in the usual finite element manner employing the developed element matrices K(e) and M(e) from Eq. (43). In the case of free damped oscillation, the linear homogenous differential equations of motion are given by

$$\begin{array}{}\text{(46)}& \mathbf{M}\ddot{\mathit{q}}\left(t\right)+\mathbf{D}\dot{\mathit{q}}\left(t\right)+\mathbf{K}\mathit{q}\left(t\right)=\mathbf{0},\end{array}$$

where MϵR^(n×n) is the system mass matrix, KϵR^(n×n) is the system stiffness matrix, DϵR^(n×n) is the system damping matrix and q(t)ϵR^(n×1) is the nodal displacement vector. The system matrices are symmetric and positively definite for finite element structures. For a free undamped system, the equation of motion is reduced to

$$\begin{array}{}\text{(47)}& \mathbf{M}\ddot{\mathit{q}}\left(t\right)+\mathbf{K}\mathit{q}\left(t\right)=\mathbf{0}\end{array}$$

By introducing the following solution approach which is given by

$$\begin{array}{}\text{(48)}& \mathit{q}\left(t\right)=\widehat{\mathit{q}}{e}^{i{\mathit{\omega }}_{\mathrm{0}}t},\ddot{\mathit{q}}\left(t\right)=\widehat{\mathit{q}}{\left(i{\mathit{\omega }}_{\mathrm{0}}\right)}^{\mathrm{2}}{e}^{i{\mathit{\omega }}_{\mathrm{0}}t}\end{array}$$

into the equation of motion (Eq. 47) the following eigenvalue problem is obtained:

$$\begin{array}{}\text{(49)}& \left({\mathbf{M}}^{-\mathrm{1}}\mathbf{K}-{\mathit{\omega }}_{\mathrm{0}k}^{\mathrm{2}}\mathbf{I}\right){\widehat{\mathit{q}}}_k=\mathbf{0},\end{array}$$

where I is a unity matrix. The condition for non-trivial solution for Eq. (49) is given by

$$\begin{array}{}\text{(50)}& p\left({\mathit{\omega }}_{\mathrm{0}k}^{\mathrm{2}}\right)=\mathit{det}\left({\mathbf{M}}^{-\mathrm{1}}\mathbf{K}-{\mathit{\omega }}_{\mathrm{0}k}^{\mathrm{2}}\mathbf{I}\right)=\mathrm{0}\end{array}$$

The n-grade characteristic polynomial $p\left({\mathit{\omega }}_{\mathrm{0}k}^{\mathrm{2}}\right)$ has n real solutions ${\mathit{\omega }}_{\mathrm{0}k},\left(k=\mathrm{1},\mathrm{\dots },n\right)$ (eigenfrequencies) and n associated eigenvectors ${\widehat{\mathit{q}}}_k$, calculated from Eq. (49). For real-life tasks the solution is usually acquired by use of eigensolver software.

5 Numerical example

The programming code for the procedure and the graphic plots described/shown below was written in MuPAD, a symbolic math toolbox in MATLAB^® (see Kusuma Chandrashekhara, 2018). The code was verified using realistic data for a wind turbine rotor blade. The blade structural data belong to a 5 MW reference wind turbine designed for offshore development (Jonkman et al., 2009). The blade is 63 m in length and is divided into 48 beam elements. The blade structural data consist of distributed mass (m_L), blade extensional stiffness (EA), flapwise stiffness (EA_zz), edgewise stiffness (EA_zz), torsional stiffness (GI_T), flapwise mass moment of inertia (Θ_y) and edgewise mass moment of inertia (Θ_z). Due to the lack of shear stiffness data in Jonkman et al. (2009), the values of (GA_sz) and (GA_sy) (used for the calculations in Table 1) – the respective edgewise and flapwise shear stiffness – are estimated to be 10 % and 20 % of extensional stiffness (EA), respectively. The values of the above-mentioned stiffness and mass moment of inertias are specified at spanwise locations along the blade pitch axis and about the principal axes of each cross section as oriented by a twist angle (γ) defined in the input data. The twist angle is incorporated using the rotational transformation of each local element stiffness and mass matrices into the global coordinate system. The results of first three (flapwise and edgewise) eigenfrequencies calculated using the Timoshenko beam model (see Kusuma Chandrashekhara, 2018), are compared with the proposed results from FAST and ADAMS in Jonkman et al. (2009). The results are as shown in Table 1.

The mode shapes and the corresponding eigenfrequencies for the first flapwise and edgewise eigenmodes as well for two torsional eigenmodes are as shown in Figs. 5–11.

In Table 2 the calculated natural frequencies for three different variants of the shear correction coefficients are shown, approximated as $\frac{G{A}_{sz}}{EA}$ or $\frac{G{A}_{sy}}{EA}$ ratios. The comparison with the frequencies calculated using the Bernoulli beam model outlines the tendency toward a stiffer structure due to the presupposed infinite shear stiffness in this case. Natural frequencies f_Bern are 0.5 %–1.0 % higher on average than f_Tim – in the (30 %, 60 %) case. The natural frequencies remain unchanged for both beam models for the purely torsional modes only. The reason for this is that the equations for torsion and bending are uncoupled (for the case of the principal axes, see Eq. 14) and remain the same in both models.

6 Conclusion and outlook

The proposed Timoshenko beam element in the 3-D description has been developed on the basis of the transfer matrix method. Both static and kinetic field matrices for the beam element are derived by applying a Runge–Kutta fourth-order numerical integration procedure on the differential equations of motion in a special way. Appropriate shape functions for the Timoshenko beam have been used to approximate the inertia forces in the motion differential equation. The beam element stiffness and mass matrices are assembled by matrix operations from the derived element field matrices. The usual finite element equations of motion for the rotor blade model are cast in the general case accounting for the structural twist angle and possible pre-bending. Therefore, in the general case the rotor blade beam model represents a polygonal approximated space curve.

For the sake of verification, the natural frequencies and associated eigenmodes are calculated using real-life rotor blade data with the incorporation of realistic twist angle data. The first two edgewise and flapwise eigenfrequencies obtained are compared with the proposed results from FAST and ADAMS software given in Jonkman et al. (2009). It can be observed that the deviation of the results of the Timoshenko beam model from FAST is comparatively lower and is in good agreement with FAST; thus, it can be stated that the approach presented regarding alternative finite element formulation works well.

One key input parameter for the Timoshenko beam model is the shear stiffness. As it was not the main goal of the present work to determine an appropriate shear correction coefficient for realistic rotor blade data, the numerical example was performed with a very rough approximation for GA_sz and GA_sy. This was done in order to simply demonstrate the performance and the differences to the Bernoulli beam model. However, if detailed data for the complex multilayer design of a rotor blade were available, a more realistic estimation of the shear stiffness could be expected. A workable method for the determination of the shear correction coefficient of a real-life rotor blade represents an important topic for further research.