Introduction to Multibody Dynamics

In the early days of the industry, wind turbine design was undertaken on the basis of quasi-static aerodynamic calculations, while the effects of structural dynamics were either ignored completely or included through the use of estimated dynamic magnification factors. From the late 1970’s research workers began to consider more reliable methods of dynamic analysis, and two basic approaches were considered: finite element representations and modal analysis.

The traditional use of standard, commercial finite element analysis software packages for solving problems of structural dynamics is challenging in the case of wind turbines. This is because of the presence of rigid body motion of one structural component, i.e. the rotor, with respect to another, i.e. the tower or another support structure type. In principle, the standard finite element method only considers structures in which the deflection occurs about an initial reference position, and for this reason the finite element models that have been developed for wind turbine in the past have been tailored to deal with this problem.

Dynamic wind turbine models commonly used as the basis of design calculations involve a modal representation of the deformation state. This approach, borrowed from the helicopter industry, has the major advantage that it offers a reliable representation of the dynamics of a wind turbine with relatively few degrees of freedom. Another important advantage is that the structural damping of flexible components can be described conveniently and realistically as modal damping. Obviously the number and type of modal degrees of freedom used to represent the dynamics of a particular wind turbine depends on the configuration and structural properties of the machine.

In principle, a wind turbine structure may be considered as a collection of a set of interconnected structural components that may undergo large rotations relative to neighbouring components, but also relative to an inertial coordinate system. The natural choice for modelling structural dynamics of said mechanical systems is the multibody dynamics approach that emerged in the late 1980’s, initially for rigid components or bodies (Nikravesh, 1988), but later also for flexible components (Shabana, 1988; Géradin and Cardona, 2001). The multibody dynamics approach is now used as an integrated part of the design process in the automotive and the aircraft industry, but it has also been used extensively in the space industry. The structural model in Bladed is based on this approach combined with a modal representation of the flexible components like the blades and the tower. This ensures a consistent formulation of the dynamic equations and facilitates the modelling of the turbines based on new structural concepts.

The multibody dynamics approach

The multibody dynamics approach used for the Bladed structural model was originally proposed by Dr. J.P. Meijaard, Nottingham University, (presently University of Twente, in the Netherlands) under commission of Garrad Hassan and Partners Ltd (Meijaard, 2005), and it was developed particularly for modelling wind turbine structural dynamics. The approach assumes a tree-like structure, which means that structural loops cannot be described outside of flexible components.

In general, the structural components may be assumed to be rigid, such as yaw and blade bearings, or flexible such as support structures, towers and blades. While rigid components are relatively easy to model, flexible components are more complicated as the motion of a material point of this type of components is generally caused by rigid body motion combined with relative motion due to the deformation. A description of the applied method for modelling flexible components is given in Modelling flexible components.

The rigid body motion of a component is described in terms of the motion of a set of local component nodes that are characteristic material points, where the motion of the component is assumed to be known. Components can only be interconnected at the nodes, and a connection is imposed by the usual finite element technique by linking the nodes of the components involved in the connection to a global structural node. Due to the assumption of the tree-like structure it is convenient to subdivide the nodes of a component into a proximal node that is linked to a node of its predecessor, and distal nodes that may be linked to nodes of successors. For example a yaw bearing has two nodes, i.e. one proximal node attached to the tower and one distal node attached to the main frame.

For all components a local Cartesian coordinate system is attached to the proximal node such that the position of origin and the orientation are defined by the position and orientation, respectively, of this node. This local coordinate system is mainly used for flexible components that are described in more details in Modelling flexible components.

The deformation state of a component is described by generalised strains that represent the degrees of freedom associated with the component. For example a yaw bearing has one generalised strain, which represent the yaw angle. The approach used also allows prescribed strains, which are particularly important in the case of stick-slip friction in bearings, where the bearing may stick if the absolute value of the angular velocity (i.e. the strain rate) approaches zero. It is noted that prescribing a strain element will reduce the effective order of the system of equilibrium equations by one.

In general, the relative motion of the distal nodes is constrained, for which reason the position and orientation of a distal node are expressed as functions of the position and orientation of the proximal node and the generalised strains. From this fundamental transformation it is straightforward to derive the corresponding transformations for the velocity and the acceleration. The constraints are conveniently expressed in terms of two constraint matrices relating to the nodal velocities and the strain rates. In general, the constraint matrices are time-dependent functions of the position and orientation of the proximal node as well as the strains.

Components may have mass or may be considered massless. The generalised inertia forces are derived from the principle of virtual work, where the inertia force density is expressed according to D’Alambert’s principle. The material velocity and acceleration are derived from a fundamental displacement hypothesis that defines the absolute displacement of all material points of the component as a function of the relative position, the nodal position and orientation and the strain. The result of this analysis shows that the inertia force can be expressed in terms of a mass matrix multiplied by the acceleration vector plus a vector of non-linear inertia forces and stresses, i.e. centrifugal and Coriolis forces.

In principle, external loads can only be applied as nodal loads or generalised stresses (equal and opposite loads applied at a generalised freedom). For distributed loads like wind and wave loads the corresponding applied nodal loads and generalised stresses are calculated according to the principle of virtual work. For a yaw bearing the applied generalised stress is simply the moment applied by the yaw actuator.

Gravity loads are conveniently considered as a distributed applied body force.

The resulting equilibrium equations of a component are derived by collecting all generalised forces acting on the component. The effect of the distal node(s) constraints is described by Lagrange’s methods in terms of internal forces that are expressed by the constraint matrices and a set of yet unknown Lagrange multipliers (Cook, Malkus and Plesha, 1989). A detailed analysis shows that the resulting component equilibrium equations and the transformation for the acceleration may be expressed in matrix form as:

\[ \begin{equation} \left[\begin{array}{ccc} \bmatrix{M}_{\mathrm{rr}}^{\mathrm{c}} & \bmatrix{M}_{\mathrm{r} \epsilon}^{\mathrm{c}} & {\bmatrix{D}_{\mathrm{r}}^{\mathrm{c}}}^T \\ {\bmatrix{M}_{\mathrm{r} \epsilon}^{\mathrm{c} }}^T & \bmatrix{M}_{\epsilon \epsilon}^{\mathrm{c}} & {\bmatrix{D}_{\epsilon}^{\mathrm{c}}}^T \\ \bmatrix{D}_{\mathrm{r}}^{\mathrm{c}} & \bmatrix{D}_{\epsilon}^{\mathrm{c}} & \bmatrix{0} \end{array}\right]\left[\begin{array}{c} \dot{\bvector{v}}^{\mathrm{c}} \\ \ddot{\bvector{\epsilon}}^{\mathrm{c}} \\ \bvector{\uplambda}^{\mathrm{c}} \end{array}\right]+\left[\begin{array}{c} \bmatrix{0} \\ \bmatrix{C}_{\epsilon \epsilon}^{\mathrm{c}} \\ \bmatrix{0} \end{array}\right] \dot{\bvector{\epsilon}}^{\mathrm{c}}+\left[\begin{array}{c} \bmatrix{0} \\ \bmatrix{K}_{\epsilon \epsilon}^{\mathrm{c}} \\ \bmatrix{0} \end{array}\right] \bvector{\epsilon}^{\mathrm{c}}=\left[\begin{array}{c} \bvector{f}_{\mathrm{a}}^{\mathrm{c}}+\bvector{f}_0^{\mathrm{c}} \\ \bvector{\upsigma}_{\mathrm{a}}^{\mathrm{c}}+\bvector{\upsigma}_0^{\mathrm{c}} \\ \bvector{a}_2^{\mathrm{c}} \end{array}\right]-\left[\begin{array}{c} \bvector{f}_{\mathrm{i}}^{\mathrm{c}} \\ \bvector{\upsigma}_{\mathrm{i}}^{\mathrm{c}} \\ \bvector{0} \end{array}\right], \label{eq:componentequilibriumequations} \end{equation} \]

where
\(\bvector{v}^{\mathrm{c}}\) is a vector of nodal velocities, \(\bvector{{\epsilon}}^{\mathrm{c}}\) is a vector of generalised and prescribed strains, \(\bvector{\uplambda}^{\mathrm{c}}\) is a vector of Lagrange multipliers corresponding to the constraints, \(\bvector{f}_{\mathrm{i}}^{\mathrm{c}}\) and \(\bvector{\upsigma}_{\mathrm{i}}^{\mathrm{c}}\) are vectors of non-linear inertia forces and stresses dual to nodal velocities and strain rates, \(\bvector{f}_{\mathrm{a}}^{\mathrm{c}}\) is a vector of applied nodal forces, \(\bvector{\upsigma}_{\mathrm{a}}^{\mathrm{c}}\) is a vector of applied generalised stresses dual to generalised strain rates, \(\bvector{f}_0^{\mathrm{c}}\) is a vector of joint reactions dual to nodal velocities, \(\bvector{\upsigma}_0^{\mathrm{c}}\) is a vector of generalised constraint stresses corresponding to prescribed strains, \(\bvector{a}_2^{\mathrm{c}}\) collects the convective terms (quadratic in nodal velocities) in the transformation for the acceleration, \(\bmatrix{M}_{\mathrm{rr}}^{\mathrm{c}}\), \(\bmatrix{M}_{\mathrm{r} \epsilon}^{\mathrm{c}}\) and \(\bmatrix{M}_{\epsilon \epsilon}^{\mathrm{c}}\) are the structural mass matrix partitions dual to nodal velocities and strain rates, \(\bmatrix{C}_{\epsilon \epsilon}^{\mathrm{c}}\) is the structural damping matrix dual to strain rates, \(\bmatrix{K}_{\epsilon \epsilon}^{\mathrm{c}}\) is the structural stiffness matrix dual to strain rates, \(\bmatrix{D}_{\mathrm{r}}^{\mathrm{c}}\) and \(\bmatrix{D}_{\epsilon}^{\mathrm{c}}\) are the constraint matrix partitions relating to nodal velocities and strain rates.

Obviously it is not possible to solve this equation due to the unknown joint reaction forces \(\bvector{f}_0^{\mathrm{c}}\) (section forces), which originates from the connection to other components. In order to solve the system it is necessary to collect all the component equilibrium equations into a global system of the complete structure, which is done using the standard finite element assembly process (Cook, Malkus and Plesha, 1989). This global system of equations has almost the same form as the component equations and is written in matrix form as

\[ \begin{equation} \left[\begin{array}{ccc} \bmatrix{M}_{\mathrm{rr}} & \bmatrix{M}_{\mathrm{r} \epsilon} & {\bmatrix{D}_{\mathrm{r}}}^T \\ {\bmatrix{M}_{\mathrm{r} \epsilon}}^T & \bmatrix{M}_{\epsilon \epsilon} & {\bmatrix{D}_{\epsilon}}^T \\ \bmatrix{D}_{\mathrm{r}} & \bmatrix{D}_{\epsilon} & \bmatrix{0} \end{array}\right]\left[\begin{array}{c} \dot{\bvector{v}} \\ \ddot{\bvector{\epsilon}} \\ \bvector{\uplambda} \end{array}\right]+\left[\begin{array}{c} \bmatrix{0} \\ \bmatrix{C}_{\epsilon \epsilon} \\ \bmatrix{0} \end{array}\right] \dot{\bvector{\epsilon}}+\left[\begin{array}{c} \bmatrix{0} \\ \bmatrix{K}_{\epsilon \epsilon} \\ \bmatrix{0} \end{array}\right] \bvector{\epsilon}=\left[\begin{array}{c} \bvector{f}_{\mathrm{a}} \\ \bvector{\upsigma}_{\mathrm{a}}+\bvector{\upsigma}_0 \\ \bvector{a}_2 \end{array}\right]-\left[\begin{array}{c} \bvector{f}_{\mathrm{i}} \\ \bvector{\upsigma}_{\mathrm{i}} \\ \bvector{0} \end{array}\right] \label{eq:componentequilibriumequationsmodified} \end{equation} \]

The main difference between component equations and the global system of equations is that the joint reaction forces do not appear in the latter as they have been cancelled out by the assembly process. Consequently the above system can be solved directly with respect to the nodal accelerations \(\dot{\bvector{v}}\), the strain accelerations \(\ddot{\bvector{\epsilon}}\), and the Lagrange multipliers \(\bvector{\uplambda}\).