Articles | Volume 7, issue 6
Research article
01 Dec 2022
Research article |  | 01 Dec 2022

A symbolic framework to obtain mid-fidelity models of flexible multibody systems with application to horizontal-axis wind turbines

Emmanuel Branlard and Jens Geisler

The article presents a symbolic framework (also called computer algebra program) that is used to obtain, in symbolic mathematical form, the linear and nonlinear equations of motion of a mid-fidelity multibody system including rigid and flexible bodies. Our approach is based on Kane's method and a nonlinear shape function representation for flexible bodies. The shape function approach does not represent the state of the art for flexible multibody dynamics but is an effective trade-off to obtain mid-fidelity models with few degrees of freedom, taking advantage of the separation of space and time. The method yields compact symbolic equations of motion with implicit account of the constraints. The general and automatic framework facilitates the creation and manipulation of models with various levels of complexity by adding or removing degrees of freedom. The symbolic treatment allows for analytical gradients and linearized equations of motion. The linear and nonlinear equations can be exported to Python code or dedicated software. There are multiple applications, such as time domain simulation, stability analyses, frequency domain analyses, advanced controller design, state observers, and digital twins. In this article, we describe the method we used to systematically generate the equations of motion of multibody systems and present the implementation of the framework using the Python package SymPy. We apply the framework to generate illustrative land-based and offshore wind turbine models. We compare our results with OpenFAST simulations and discuss the advantages and limitations of the method. The Python implementation is provided as an open-source project.

1 Introduction

The next generation of wind turbine digital technologies and control systems require versatile aero-servo-hydro-elastic models, with various levels of fidelity, suitable for a wide range of applications. Such applications include time domain simulations, linearization (for controller design and tuning, or frequency domain analyses), analytical gradients (for optimization procedures), and generation of dedicated, high-performance or embedded code (for stand-alone simulations, state observers or digital twins). Current models are implemented for a specific purpose and are usually based on an heuristic structure. Aeroelastic tools, such as Flex (Øye1983; Branlard2019) or ElastoDyn (Jonkman et al.2021), rely on an assumed chain of connections between bodies, a given set of degrees of freedom, and predefined orientations of shape functions. It is not straightforward to extract reduced-order models from these tools or extend the models to additional degrees of freedom.

Tools with linearization capabilities, such as HAWCStab2 (Sønderby and Hansen2014) or OpenFAST (Jonkman et al.2021), are dedicated to horizontal-axis wind turbines, and the evaluation of the gradients are limited to hard-coded analytical expressions or numerical finite differences. Small implementation changes often require extensive redevelopment, and the range of applications of the tools remains limited (Simani2015). The linear models generated by these tools are numerical models that are evaluated for a given set of numerical input parameters. It is therefore difficult to obtain gradients of the linear models as function of the input parameters, information which is becoming increasingly important in optimization frameworks and controls co-design approaches (Jonkman et al.2022).

To address these issues, we propose a symbolic framework (also called a computer algebra program) for the automatic derivation, processing, and parameterization of models with granularity in the level of fidelity. Our approach is based on Kane's method (Kane and Wang1965) and a nonlinear shape function representation of flexible bodies (Shabana2013) described using a standard input data (SID) format (Wallrapp1994; Schwertassek and Wallrapp1999). The method yields compact symbolic equations of motion with implicit account of the constraints. Similar approaches have been presented in the literature: Kurz and Eberhard (2009), Merz (2018), Lemmer (2018), and Branlard (2019). Our framework differs in the fact that all equations are processed at a symbolic level, and therefore the model can be used in its nonlinear or linearized form. The linear models are obtained using analytical differentiation. They can be evaluated for various sets of input parameters directly and therefore be used in optimization frameworks or control co-design approaches. We implemented an open-source version in Python using SymPy (SymPy2021), leveraging its mechanical toolbox. Alternative symbolic frameworks found in the literature are usually limited to rigid bodies (Verlinden et al.2005; Kurz and Eberhard2009; Gede et al.2013; Docquier et al.2013) or are closed-source or using proprietary software (Reckdahl and Mitiguy1996; Kurtz et al.2010; MotionGenesis2016; Lemmer2018).

Kane's method and the nonlinear shape function approach presented in this article do not represent the state of the art of multibody dynamics with flexible bodies. The geometrically exact beam theory (Simo1985; Jelenić and Crisfield1999; Géradin and Cardona2001; Bauchau2011) is more precise than the shape function approach because it represents the beam kinematics exactly. Linearization of the geometrical exact beam theory equations is possible and also more precise than the shape function approach, but it leads to larger and more involved expressions. Similarly, multipurpose multibody software exists (Lange et al.2007), such as ANSYS (ANSYS2022), SIMPACK (SIMPACK2022), or MBDyn (MBDyn2022). These more advanced approaches target different applications than those envisioned in this study: they are suitable for numerical simulations, but they cannot provide symbolic mid-fidelity models in compact form.

In Sect. 2, we present the formalism used to derive the equations of motion in a systematic and unified way for flexible and rigid bodies. In Sect. 3, we give an overview of how the symbolic calculation framework was implemented. Example of applications relevant to wind energy are given in Sect. 4. Discussions and conclusions follow.

2 Method to obtain the equations of motion

In this section, we present the formalism used to set up the equations of motion.

2.1 System definition and kinematics

We consider a system of nb bodies, rigid or flexible, connected by a set of joints. For simplicity, we assume that no kinematic loops are present in the system and the masses of the bodies are constant. An inertial frame is defined to express the positions, velocities, and accelerations of the bodies. We adopt a minimal set of generalized coordinates, q, of dimension nq, to describe the kinematics of the bodies: joint coordinates describing the joints displacements and Rayleigh–Ritz coordinates for the amplitudes of the shape functions of the flexible bodies (see, e.g., Branlard2019). The choice of coordinates is left to the user, but it is assumed to form a minimal set. We will provide illustrative examples in Sect. 4.

At a given time, the positions, orientations, velocities, and accelerations of all the points of the structure are entirely determined by the knowledge of q,q˙, and q¨, where (˙) represents the time derivative. For a given body i and a point P belonging to the body, the position, velocity, and acceleration of the point are given by (see, e.g., Shabana2013)


where ri, vi, and ai are the position, velocity, and acceleration of the origin of the body, respectively; sP0 is the initial (undeformed) position vector of point P with respect to the body origin; the subscript P is used for the deformed position of the point and P0 for the undeformed position; uP is the elastic displacement of the point (equal to 0 for rigid bodies); ωi is the rotational velocity of the body with respect to the inertial frame; and (˙) and (˙)i refer to time derivatives in the inertial and body frame respectively. Throughout the article, we use bold symbols for vectors and matrices and uppercase symbols for most matrices. The elastic displacement is obtained as a superposition of elastic deformations (see Sect. 2.4). We define the transformation matrix Ri that transforms coordinates from the body frame to the inertial frame, and by definition [ω̃i]=R˙iRiT, where [̃] represents the skew symmetric matrix, and the exponent T denotes the matrix transpose. We assume that vectors are represented as column vectors to conveniently introduce matrix-vector multiplications. We use the notation “” to indicate the dot product between two vectors (irrespective of their column or row representation).

2.2 Introduction to Kane's method

Kane's method (Kane and Wang1965) is a powerful and systematic way to obtain the equations of motion of a system. The procedure leads to nq coupled equations of motion:

(4) f r + f r = 0 , r = 1 n q ,

where fr is associated with inertial loads and fr is associated with external loads, and these components are obtained for all generalized coordinates. The components are obtained as a superposition of contributions from each body:

(5) f r = i = 1 n b f ri , f r = i = 1 n b f ri .

The terms fri and fri can be obtained for each body individually and assembled at the end to form the final system of equations. We will present in Sects. 2.3 and 2.4 how these terms are defined for rigid bodies and flexible bodies, respectively.

2.3 Rigid bodies

We assume that body i is a rigid body and proceed to define the terms fri and fri. The inertial force, fi, and inertial torque, τi, acting on the body are

(6) f i = - m i a G , i , τ i = - I G , i ω ˙ i - ω i × ( I G , i ω i ) ,

where mi is the mass of the body, aG,i is the acceleration of its center of mass with respect to the inertial frame, and IG,i is the inertial tensor of the body expressed at its center of mass. Equation (6) is a vectorial relationship; it may therefore be evaluated in any coordinate system. The component fri is defined as

(7) f ri = J v , ri f i + J ω , ri τ i ,


(8) J v , ri = v G , i q ˙ r , J ω , ri = ω i q ˙ r ,

where vG,i is the velocity of the body mass center with respect to the inertial frame. The partial velocities, or Jacobians, Jv and Jω, are key variables of Kane's method. They project the physical coordinates into the generalized coordinates (q), inherently accounting for the kinematic constraints between bodies. In numerical implementations, the Jacobians are typically stored in matricial forms, referred to as “velocity transformation matrices”. The terms fri can equivalently be obtained using the partial velocity of any body point (e.g., the origin) by carefully transferring the inertial loads to the chosen point.

The external forces and torques acting on the body are combined into an equivalent force and torque acting at the center of mass, written as fi and τi. The component fri is then given by

(9) f ri = J v , ri f i + J ω , ri τ i .

Equivalently, the contributions from each individual force, fi,j, acting on a point Pj of the body i, and each torque, τi,k, can be summed using the appropriate partial velocity to obtain fri:

(10) f ri = j v P j q ˙ r f i , j + k J ω , ri τ i , k ,

where vPj is the velocity of the point j with respect to the inertial frame. Equations (7) and (9) are inserted into Eq. (5) to obtain the final equations of motion.

2.4 Flexible bodies

We assume that body i is a flexible body and proceed to define the terms fri and fri. The dynamics of a flexible body are described in standards textbooks such as Shabana (2013) or Schwertassek and Wallrapp (1999). Unlike rigid bodies, the equations for flexible bodies are typically expressed with respect to a reference point different from the center of mass. We will call this point the origin and write it Oi. The elastic displacement field of the body is written as u. It defines the displacement of any point of the body with respect to its undeformed position. Using the zeroth-order1 Rayleigh–Ritz approximation, the displacement field at a given point, P, is given by the sum of shape function contributions:

(11) u ( P ) = j = 1 n e , i Φ i j ( P ) q e , i j ( t ) ,

where Φij are the shape functions (displacement fields) of body i, and qe,ij is the subset of q consisting of the elastic coordinates of body i, of size ne,i. The principles of the shape function approach applied to beams are given in Appendix B. The shape functions are more easily represented in the body coordinate system. Vectors and matrices that are explicitly written in the body frame will be written with primes. The equations of motion of the flexible bodies are (Wallrapp1994)

(12) M xx M x θ M xe M θ θ M θ e sym. M ee i a i ω ˙ i q ¨ e , i + k ω , x k ω , θ k ω , e i + 0 0 k e i = f x f θ f e i ,

where the x, θ, and e subscripts respectively indicate the translation, rotation, and elastic components; M is the mass matrix of dimension 6+ne,i made of the block matrices Mxx,,Mee; ai and ω˙i are the linear and angular acceleration of the body (origin) with respect to the inertial frame; kω represents the centrifugal, gyration, and Coriolis loads, also called quadratic velocity loads; ke represents the elastic strain loads, which may contain geometric stiffening effects; and f represents the external forces, torques, and elastic generalized forces. The different components of M, kω, ke, and f are given in Appendix A. These terms depend on q, q˙, and Φi. The inertial force, torque, and elastic loads are


The external and elastic loads are


The components of fri and fri, for r=1…nq, are then defined as



(21) J v , ri = v O , i q ˙ r , J ω , ri = ω i q ˙ r , J e , ri = q e , i q r ,

where vO,i is the velocity of the body with respect to the inertial frame. The term Je,ri consists of 0 and 1 because qe,i is a subset of q. Equations (19) and (20), once evaluated for body i, are inserted into Eq. (5) to obtain the final equations of motion.

2.5 Nonlinear and linear equations of motion

The nq equations of motion given in Eq. (4) are gathered into a vertical vector e. They are recast into the form

(22) e ( q , q ˙ , q ¨ , u , t ) = f + f = F ( q , q ˙ , u , t ) - M ( q ) q ¨ = 0


(23) M ( q ) q ¨ = F ( q , q ˙ , u , t ) ,

where M=-eq¨ is the system mass matrix and F is the forcing term vector – that is, the remainder terms of the equation (F=e+Mq¨). In this section, we introduce the vector u to represent the time-dependent inputs that are involved in the determination of the external loads (not to be confused with the displacement field introduced in Eq. 11). Both sides of the equations are also dependent on some parameters, but this dependency is omitted to shorten notations. The stiffness and damping matrices may be obtained by computing the Jacobian of the equations of motion with respect to q and q˙, respectively. The nonlinear equation given in Eq. (23) is easily integrated numerically, for instance by recasting the system into a first-order system or by using a dedicated second-order system time integrator.

In various applications, a linear time-invariant approximation of the system is desired. Such approximation is obtained at an operating point, denoted with the subscript 0, which is a solution of the nonlinear equations of motion, namely

(24) e ( q 0 , q ˙ 0 , q ¨ 0 , u 0 , t ) = 0 .

The linearized equations about this operating point are obtained using a Taylor series expansion:

(25) M 0 ( q 0 ) δ q ¨ + C 0 ( q 0 , q ˙ 0 , u 0 ) δ q ˙ + K 0 ( q 0 , q ˙ 0 , q ¨ 0 , u 0 ) δ q = Q 0 ( q 0 , q ˙ 0 , u 0 ) δ u ,


(26) M 0 = - e q ¨ 0 , C 0 = - e q ˙ 0 , K 0 = - e q 0 , Q 0 = e u 0 ,

where M0, C0, and K0 are the linear mass, damping, and stiffness matrices, respectively; Q0δu is the linear forcing vector (Q0 is the input matrix); δ indicates a small perturbation of the quantities; and |0 indicates that the expressions are evaluated at the operating point. In practical applications, linearization is done at an operating point where the acceleration is zero (q¨0=0) and most velocities are also zero. Examples of applications of the linear equations of motion are controller design, frequency domain analyses, and stability analyses. The symbolic system matrices also allow for the easy formulation of linear parameter-varying models used in many advanced control applications.

3 Implementation into a symbolic framework

In this section, we discuss the Python open-source symbolic calculation framework that we implemented according to the equations given in Sect. 2. A Maxima implementation from the same authors is also available (Geisler2021).

The Python library YAMS (Yet Another Multibody Solver) started as a numerical tool published in previous work (Branlard2019). The library is now supplemented with a symbolic module so that both numerical and symbolic calculations can be achieved. The new implementation uses the Python symbolic calculation package SymPy (SymPy2021). We leveraged the features present in the subpackage “mechanics”, which contains all the tools necessary to compute kinematics: the definition of frames and points and the determination of positions, velocities, and accelerations. The subpackage also contains an implementation of Kane's equations for rigid bodies (i.e., Sect. 2.3). We were also inspired by the package PyDy (Gede et al.2013), which is a convenient tool to export the equations of motion to executable code and directly visualize the bodies in 3D. The core of our work consisted of implementing a class to define flexible bodies (FlexibleBody) and the corresponding Kane method for this class (Sect. 2.4).

For the FlexibleBody class, we followed the formalism of Wallrapp (1994) and implemented Taylor expansions for all the terms defined in Appendix A, allowing the symbolic computation with Taylor expansions to any order. In practice, a zeroth- or first-order expansion is used. The use of Taylor expansions is presented in Appendix D3. The different Taylor coefficients may be kept as symbolic terms or replaced early on by numerical values provided by a SID, for instance.

We structured the code into three layers.

  1. The low-level layer integrates seamlessly with SymPy and PyDy by using the FlexibleBody class we provide. It is the layer that offers the highest level of granularity and control for the user, since arbitrary systems with various kinematic constraints can be implemented, at the cost of requiring more expertise.

  2. The second level automates the calculation of the kinematics by introducing simple connections between rigid and flexible bodies. The connections may be rigid, with constant offsets and rotations, or dynamic. A connection from a flexible body to another body is assumed to occur at one extremity of the flexible body. Some knowledge of SymPy mechanics is still required to use this layer.

  3. The third level consists of template models such as generic land-based or offshore wind turbine models. Degrees of freedom are easily turned on and off for these conceptual models depending on the level of fidelity asked by the user, and generic external forces can be implemented or declared as external inputs.

Figure 1Typical workflow for the usage of the symbolic framework, going from numerical inputs and a conceptual model to numerical packages that can be used for various applications.


The overall workflow for typical usage of the symbolic framework is illustrated in Fig. 1. The symbolic framework takes as input a conceptual model of the structure, which is assembled using one of the three layers previously described. The nonlinear and linear equations of motion can be exported to LaTeX and Python-ready scripts for various applications (see Sect. 5.1). Using the third layer, as little as three lines of code are required by the user to perform the full step from derivation of the equations, optional linearization, and exportation. To obtain numerical results from the exported Python code, the user needs to provide the arrays with the degrees of freedom values q and q˙, their initial conditions, a dictionary with inputs (u) that are functions of time, and a dictionary of parameters (p) containing all the numerical constants such as mass, acceleration of gravity, and geometric parameters. We implemented various preprocessing tools in YAMS to facilitate the calculation of numerical parameters, typically from a set of OpenFAST input files or by using structural parameters defined by the users. YAMS contains tools to compute the flexible bodies parameters (mass matrix, stiffness matrix, shape integrals) using integrals over the shape functions or using a finite-element beam formulation. YAMS also contains tools to compute the rigid body inertia of different components of a wind turbine or the full system. Postprocessing tools are also included to readily time-integrate the generated model using numerical values (including initial values).

The source code of YAMS is available on GitHub as a subpackage of the Wind Energy Library, WELIB (Branlard2022a). The repository contains tests and working examples, including the ones presented in Sect. 4.

4 Wind energy applications

4.1 Approach

In this section, we present different wind energy applications of the symbolic framework. We focus on models with at least one flexible body because the rigid body formulation of SymPy has been well verified (Gede et al.2013). For each example, the equations of motion are given and their results are compared with OpenFAST (Jonkman et al.2021) simulations. This is readily achieved because our framework can export the equations of motion to Python functions, load input files from an OpenFAST model, and integrate the generated equations using the same conditions as defined in the OpenFAST input files. In this article, we do not focus on the modeling of the external loads, but we include them in the equations of motion. It is the responsibility of the user to define these functions, for instance through aero- or hydro-force models. For the verification results presented in this section, we only include the gravitational and inertial loading. In all examples, the National Renewable Energy Laboratory (NREL) 5 MW reference wind turbine (Jonkman et al.2009) is used. The examples below are provided on the GitHub repository where the YAMS package is provided (Branlard2022a).

4.2 Notations

We adopt a system of notations where the first letter of a body is used to identify the parameters of that body. As an example, the tower is represented with the letter T, and the following body parameters are defined: T, origin; MT, mass; LT, length; (Jx,T,Jy,T,Jz,T), diagonal coefficients of the inertia tensor about the center of gravity and in body coordinates; rTG, vector from body origin to body center of mass of coordinates (xTG,yTG,zTG) in body coordinates. We also define θt, the nacelle tilt angle about the y axis; g, the acceleration of gravity along z; and O, the origin of the global coordinate system.

4.3 Rotating blade with centrifugal stiffening

We begin with the study of a flexible blade of length LB=R, rotating at the constant rotational speed Ω. We use this test case to familiarize the reader with the key concepts of the shape function approach given in Appendix B. A sketch of the system is given in Fig. 2. We model the blade using a single shape function, assumed to be directed along the x axis (in the “flapwise” direction): Φ1=Φ1,xex=Φex, where ex is the unit vector in the x direction, and Φ (without subscripts) is used to shorten notations. The undeflected blade is directed along the radial coordinate r and rotates around the x axis. We assume that the shape function is known, denoted Φ(r). It can be computed as the first flapwise mode of the blade using tools provided in YAMS. The expression Φ(r)=(r/R)3 is a simple approximation that can be used for hand calculations. The displacement at a given radial position and a given time is given by Eq. (11): u(r,t)=Φ1(r)qe,1(t)=Φ(r)q(t)ex, where qe,1 is the generalized coordinate associated with Φ1, and we use q without subscripts to shorten notations. The aerodynamic force per length in the flapwise direction is denoted px(r). The generalized mass and stiffness are computed based on the mass per length (m) and flapwise bending stiffness (EIy) of the blade, according to Eq. (B1):


Figure 2Sketch of a rotating blade with the restoring centrifugal force. Points are indicated in green, degrees of freedom in blue, and loads in orange.


The generalized force is obtained from Eq. (B3):

(29) f e = 0 R p x ( r , t ) Φ ( r ) d r .

The important consideration for this model is the axial load, N. The main axial load at a radial station r comes from the centrifugal force acting on all the points outboard of the current station:

(30) N ( r ) = r R m ( r ) Ω 2 r d r .

The geometric stiffness contribution of the axial load is obtained from Eq. (B5) as

(31) K g ( Ω ) = 0 R N ( r ) d Φ d r 2 d r = Ω 2 0 R r R m ( r ) r d r d Φ d r 2 d r .

The geometric stiffness, Kg, is positive and increases with the square of the rotational speed. This restoring effect is referred to as “centrifugal stiffening”. In this example, the beam rotates with respect to a fixed support, the influence of gravity is omitted, and no force other than the centrifugal force is assumed in the radial direction (the Coriolis force contribution to the radial force is assumed to be negligible for simplicity). Therefore, the only geometric stiffness comes from the centrifugal force. For a wind turbine blade mounted on a flexible support and under the influence of gravity, the different geometric stiffening terms presented in Appendix C should be used. Adding the elastic and geometric stiffness, the natural frequency of the blade increases with the rotational speed as follows:

(32) ω 0 ( Ω ) = ( K e + K g ( Ω ) ) M e = ω 0 2 ( 0 ) + K g ( Ω ) M e = ω 0 2 ( 0 ) + k Ω Ω 2 ,

where kΩ is referred to as the “rise factor” or “Southwell coefficient”, and in our approximation, it is found to be constant: kΩ=Kg(Ω)/Me/Ω2. The coefficient provides the variation of the blade frequency with rotational speed, which is something that is observed on a Campbell diagram when performing stability analyses. In general, the mode shapes of the blade will also change as a function of the rotational speed, and different shape functions should preferably be used for simulations at different rotational speeds. The effect is fairly limited, and most OpenFAST practitioners only use one shape function corresponding to the value at rated rotational speed. Similarly, the Southwell coefficient is a function of the rotational speed, but the variation is negligible as long as the rotational speed is small compared to the natural frequency (e.g., (Ω/ω)25; see Bielawa2006), which is the case for wind energy applications.

The treatment for a shape function purely in the blade edgewise direction is similar, using Φ22eθ, where eθ is the unit vector in the edgewise direction. In this case, the centrifugal force also has a component in the tangential direction, pθ,centri(r)=-Ω2uθ(r)dm(r), with uθ2q2. This leads to a generalized force equal to 0Lpθ,centriΦ2dr=-Ω2Meq2, or, equivalently, to a stiffness term: Kω=-Ω2Me. It can be verified that this generalized force corresponds to the contribution Oe,11ωx2, from kω,e, given in Eq. (A10). For an edgewise mode, the frequency therefore evolves as

(33) ω 0 ( Ω ) = ( K e + K g ( Ω ) + K ω ( Ω ) ) M e = ω 0 2 ( 0 ) + ( k Ω - 1 ) Ω 2 ,

with kΩ=Kg(Ω)/Me/Ω2 and with Kg computed using Eq. (31).

We apply the method to the NREL 5 MW wind turbine using the blade properties and shape functions provided in the ElastoDyn input file. We order the degrees of freedom as first flap, first edge, and second flap, assuming no coupling between the shape functions, so that each can be treated individually using the results from this section. The diagonal coefficients of the mass matrix are diag(Me)=[9.5×103,1.5×104,5.7×103], and for the stiffness matrix they are diag(Ke)=[1.7×104,6.7×104,8.7×104], computed according to Eqs. (27) and (28). The coefficients kΩ of each degree of freedom are obtained as kΩ=[1.7,1.4,5.5]. We compare the frequencies obtained with the present method against OpenFAST linearization results in Fig. 3. The simulations were run in vacuum (no gravity, no aerodynamics) and with a cone angle of 0. Strong agreement is found for the evolution of the different frequencies with the rotational speed. The stiffening is less pronounced for edgewise modes as a result of the softening introduced by Kω.

Figure 3Variation of the natural frequencies of the NREL 5 MW turbine blade with rotational speed. Results from YAMS and OpenFAST, with mean relative error, ϵ, are reported on the figure.


This section focused on the analysis of individual shape functions, expressed purely in one direction. In the general case, multiple shape functions are present, and couplings might exist between them (due to the structural twist or nonorthogonality of the shape functions, or if the shape functions have components in multiple directions). For instance, if the blade is represented by two shape functions defined in both the x and y direction, we have Φ1=Φ1,xex+Φ1,yey, Φ2=Φ2,xex+Φ2,yey, and u=Φ1q1+Φ2q2. The developments of appendices A and B should be used for general cases.

Figure 4Model of a land-based or fixed-bottom wind turbine using 1 to 3 degrees of freedom (fore–aft and side–side flexibility of the support structure, as well as shaft rotation). Points are indicated in green, degrees of freedom in blue, and loads in orange.


4.4 Two-degrees-of-freedom model of a land-based or fixed-bottom turbine

We consider a system of three bodies: tower (or support structure), nacelle, and rotor. The system represents a land-based wind turbine or a fixed-bottom offshore wind turbine. A sketch of the system is given in Fig. 4. The nacelle and rotor blades are rigid bodies, whereas the tower is flexible and represented by one shape function2 in the fore–aft direction, denoted Φ11ex. For hand calculations and as a first approximation, the first mode shape of a massless beam with a top mass may be used: Φ1(z)=1-cos(zπ/L/2). Increased accuracy is obtained when the shape function matches the actual first tower fore–aft bending mode, accounting for the effect of the rotor–nacelle mass and inertia. The degrees of freedom are q=(q,ψ), where q is the generalized (elastic) coordinates in the fore–aft direction and ψ is the azimuthal position. The slope of the tower shape function at the tower top is a key coupling parameter of the model, denoted νy. When the tower deflects 1 m in the x direction, the nacelle rotates by an angle νy. The method assumes that the tower-top point remains along the x axis, neglecting the so-called nonlinear geometric effect. However, nonlinear geometric effects can be included using geometric stiffening corrections (see Appendix C or Branlard2019). The aerodynamic thrust and torque are denoted fa and τa, respectively, and act at the rotor center (point R). The low-speed shaft generator torque is written as τg. The distributed loads on the tower, px (from aerodynamics and hydrodynamics), are projected against the shape function to obtain the generalized forces fe=0LTpx(z,t)Φ1(z)dz. The moments of inertia of the rotor in its coordinates are (Jx,R,J,R,J,R). We note that Me,Ke, and De are the generalized mass, stiffness, and damping, respectively, associated with a given shape function Me=0LTm(z)Φ12(z)dz, Ke=0LTEI(z)d2Φ1dz2(z)2dz, De=2ζMeωe, where m(z) and EI(z) are the mass per length and bending stiffness of the tower, respectively, and ωe and ζ are the frequency and damping ratio, respectively, associated with the shape function (assuming the shape function approximates a mode shape). The geometric softening of the tower due to the tower-top mass (Kgt) and its own weight (Kgw) is obtained using Eq. (B5), as Kg=Kgt+Kgw, with


The tower is assumed to be fixed and under no significant vertical external loads, and therefore the only geometric stiffness comes from the gravitational force. For a tower mounted on a moving support (fixed-bottom foundation or floater), additional geometric stiffening terms would be present (see Appendix C). The shape function frequency is obtained as

(36) ω e = ( K e + K g ) / M e .

The application of the symbolic framework leads to the following equations of motion (rearranged for interpretability):

(37) M q 0 0 J x,R q ¨ ψ ¨ = f q τ a - τ g ,





Details on the derivations are given in Appendix E1. The mass matrix consists of three main contributions: Eq. (38) represents the elastic mass and the rotor nacelle assembly (RNA) mass, Eq. (39) is the generalized rotational inertia of the RNA, and Eq. (40) is the inertial coupling between the tower bending and the rotation of the nacelle. The forcing terms are identified as follows: Eq. (41) consists of the elastic load resulting from the external forces on the tower, the elastic and geometric stiffness loads, and the damping load on the tower; Eq. (42) is the gravitational load from the RNA, which will contribute to the stiffness of the system; Eq. (43) is the centrifugal force of the RNA (Mω2r with ω=νyq˙); Eq. (44) is the generalized torque from the aerodynamic thrust; and Eq. (45) is the thrust contribution acting directly along the direction of the shape function degree of freedom (along x). The RNA center of mass plays an important part in the equations (see the terms (MNxNG+MRxNR) and (MNzNG+MRzNR)).

The equations of motion given in Eq. (37) can be used to perform time domain simulations of a wind turbine. It is noted that the 2 degrees of freedom are only coupled by the aerodynamic loads. The nonlinear model was used in previous work for time domain simulations, and its linear version was used for state estimations (Branlard et al.2020a, b). In this section, we apply the linearized form to compute the natural frequency of the turbine tower fore–aft mode. The linearized stiffness is obtained by taking the gradient of the forcing with respect to q and using a small angle approximation for νy to the second order:

(46) K q , lin = ( K e + K g ) - ν y 2 g ( M N z NG + M R z NR - f a q cos θ t ) + ν y f a sin θ t .

For the NREL 5 MW reference turbine (Jonkman et al.2009), the different numerical values are g=9.807m s−2, θt=5, xNR=-5.0m, zNR=2.4m, LT=87.6m, zNG=1.75m, xNG=1.9m, MR=1.1×105kg, Jx,R=3.86×107kg m2, J,R=1.92×107kg m2, MN=2.4×105kg, Jy,N=1.01×106kg m2, and MRNA=3.5×105kg. The first fore–aft shape function of the NREL 5 MW turbine tower and its derivatives are


with z=z/L, a2=0.7004, a3=2.1963, a4=-5.6202, a5=6.2275, and a6=-2.504. The material properties and the shape function are illustrated in Fig. 5. The scaling of the shape functions given in Eq. (47) is important to obtain the correct numerical values for the flexible tower, namely νy=0.0185, Me=5.4×104, Ke=1.91×106, Kg=-5.2×104-1.0×104=-6.20×104, ωe=(Ke+Kg)/Me=5.85rad s−1. These numerical values, with q=0, lead to Mq=4.375×105 and Kq=1.849×109. The first fore–aft mode of the wind turbine has a natural frequency of f=Kq/Mq=0.3272Hz. This value was compared with results obtained using OpenFAST linearization. Both methods are in strong agreement, with differences only arising at the fifth decimal place.

4.5 Three-degrees-of-freedom model of a land-based or fixed-bottom turbine

We consider the same system as the one presented in Sect. 4.4, but the tower is now represented by one shape function in both the fore–aft and side–side directions, Φ11ex and Φ22ey. The degrees of freedom are q=(q1,q2,ψ), where q1 and q2 are the generalized (elastic) coordinates in the fore–aft and side–side directions, respectively, and ψ is the rotor azimuth. A sketch of the system is given in Fig. 4.

Figure 5Properties of the NREL 5 MW turbine tower: mass per length (m), bending stiffness (EI), shape function displacement (Φ), slope (dΦ/dz), and curvature (d2Φ/dz2).


The slopes of the shape functions at the tower top are key coupling parameters of the model, denoted νx and νy. The aerodynamic thrust and torque are denoted fa and τa, acting at point R. The distributed loads on the tower, px and py (from aerodynamics and hydrodynamics), are projected against the shape functions to obtain the generalized forces fe1=Φ1pxdz and fe2=Φ2pydz. The moments of inertia of the rotor in its coordinates are (Jx,R,J,R,J,R). We note that Me, Ke, and De are the generalized mass, stiffness, and damping, respectively, associated with a given shape function (e.g., Me11=Φ12m(z)dz, where m is the mass per length of the tower). The application of the symbolic framework leads to the equations of motion given in Appendix E2. To simplify the equations and limit their length when printing them in this article, we have applied a first-order small-angle approximation for θt and a second-order approximation for νx and νy. It is observed from Eq. (E14) that a first-order approximation for νy would have removed the influence of the rotor and nacelle y inertia on the generalized mass associated with the tower fore–aft bending.

We performed a time simulation of the model using both our symbolic framework YAMS and OpenFAST. The time integration in YAMS currently relies on tools provided in the SciPy package, which implements several time integrators. A sufficient level of accuracy was obtained using a fourth-order Runge–Kutta method, which is the default method. Kane's method, which uses a minimal set of coordinates, tends to lead to stiff systems, and it is possible that implicit integrators may be needed for other systems. We compare the time series obtained using our generated functions with results from the equivalent OpenFAST simulation in Fig. 6. In this simulation, the tower top is initially displaced by 1m in the x and y directions, and the rotational speed is 5rpm. We report the mean relative error, ϵ, and the coefficient of determination, R2, on the figure. We observe that our model is in strong agreement with the OpenFAST simulation. The differences in the second tower degree of freedom are attributed to (1) the handling of the small-angle approximation, which is different in OpenFAST (using the closest orthonormal matrix; Jonkman2009) and in our formulation (two successive rotations, linearized), and (2) the nonlinear geometric corrections that are implemented in OpenFAST, which we have omitted here by only selecting shape function expansion to the zeroth order (see Sect. 5.2). The variation in azimuthal speed, resulting from the coupling between the gyroscopic loads and the tower bending, is captured well.

Figure 6Free decay results for the land-based/fixed-bottom model using both the symbolic framework (YAMS) and OpenFAST. From top to bottom: tower fore–aft bending, tower side–side bending, and shaft rotational speed.


4.6 Three-degrees-of-freedom model of a floating wind turbine

In this example, we demonstrate the applicability of the method for a floating wind turbine. We model the turbine using three bodies: rigid floater, flexible tower, and rigid RNA (labeled “N”). The degrees of freedom selected are q=(x,ϕ,qT), where x is the floater surge, ϕ is the floater pitch, and qT is the coordinate associated with a selected fore–aft shape function. A sketch of the model is given in Fig. 7. The notations are similar to the ones presented in Sect. 4.5. Lumped hydrodynamic loads at the floater center of mass are now added. The model can also be used for a combined tower and floater that is flexible, simply by setting the mass of the floater to zero and including the hydrodynamic loading into the loading px. The equations of motion are given in Appendix E3. The equations were simplified using a first-order small-angle approximation of θt and ϕy and a second-order approximation for νy.

Figure 7Model of a floating wind turbine using 3 degrees of freedom. Points are indicated in green, degrees of freedom in blue, and loads in orange.


We performed a numerical simulation of the model generated by YAMS and compared it with OpenFAST for a case with gravitational loads only, starting with x=0m, ϕ=2, and qT=1m. The results are presented in Fig. 8. We observe again that the results from the two models correlate to a high degree.

Figure 8Free-decay results for the floating wind turbine model using YAMS and OpenFAST. From top to bottom: surge, pitch, and tower fore–aft bending.


We also compared the linearized version of both models. The symbolic framework can generate the linearized mass, stiffness, and damping matrices, as described in Sect. 2.5. The matrices are then combined into a state matrix and compared with the state matrices written by the OpenFAST linearization feature. The eigenvalue analysis of the YAMS state matrix returned a pitch and fore–aft frequencies of 0.099 and 0.799Hz, respectively, whereas OpenFAST returned 0.095 and 0.795Hz. The 4 % error in the pitch frequency appears reasonable in view of the approximations used.

5 Discussions

5.1 Applications and advantages of the method

The implementation of the symbolic YAMS library was originally motivated by the need to obtain a simple linearized model of a floating wind turbine for frequency domain simulations. There are multiple potential applications of the framework.

  • The generated equations can be used in time domain simulation tools. The equations can be readily exported to different programming languages (C, FORTRAN, or Python) providing computationally efficient tools, particularly because the method generates compact and minimal equations. This is in contrast to most other multibody codes, in which many terms are calculated as matrix equations and through successive function calls. Further, the symbolic framework allows us to generate optimized code, in which common terms and factors are computed once and stored in temporary variables for reuse in the different expressions. In our examples, time domain simulations were observed to be 2 orders of magnitude faster when using the automatically generated code in Python compared to OpenFAST simulations that rely on a compiled language. Using such a framework can be considered in the future to replace the existing ElastoDyn module of OpenFAST. It can also be applied to unusual configurations such as multirotor or vertical-axis turbine concepts. Dedicated code can be generated for specific applications for increased performance. For instance, implicit integrators with iterative Newton–Raphson-like solvers benefit from the possibility of generating exact and efficient Jacobians along with the equations of motion.

  • The generation of linearized models has a wide range of applications, such as linear time domain simulations, controller design and tuning, frequency domain analyses, stability analyses, state observers, or digital twins. The symbolic approach is severalfold faster than alternative approaches because it can be evaluated for all operating points at once, whereas other methods (e.g., OpenFAST, HAWCStab2) require multiple linearization calls.

  • Analytical linearization with respect to parameters is directly obtained using our tool, which can be used for sensitivity analyses, parameter studies, optimizations, integrated design approaches, and controls co-design (e.g., using methods such as linear-matrix-inequality-based designs; Pöschke et al.2020). Nonanalytical approaches require numerous linearizations and evaluations at various operating points (Jonkman et al.2022).

  • In addition to the nonlinear or linear equations of motion in minimal coordinates, the equations for the constraint forces or any auxiliary kinematic variable can also be generated efficiently by inserting unknown virtual displacements in the equations (see Appendix D5 for an alternative approach). The position of all bodies in local or global coordinates can be recovered from the minimal coordinates and, in combination with the flexible code generation, be used to output data (e.g., for 3D animations of the turbine).

  • Analytical gradients of the equations can be computed and used in optimizations, nonlinear model predictive control, or moving horizon estimation. External loads that cannot be expressed analytically can be defined as generic functions of the structural degrees of freedom, inputs, and parameters. After the code generation, the user can link a numerical implementation of the function and its numerical gradients to be able to use a mix of analytical and numerical gradients.

  • Another advantage of the presented method is the possibility to quickly generate models with different levels of detail, ensuring consistency between the different levels of fidelity. This is in contrast to other more heuristic modeling approaches in which parameters often have to be retuned for each added degree of freedom.

  • The method provides useful insights and can be used as an educational tool: simple models of a system with few degrees of freedom can readily be obtained, studied, and compared to hand-based calculation.

5.2 Advanced consideration

Section 2 addressed the systematic derivation of the equations of motion for an assembly of rigid or flexible bodies. Some advanced aspects of the method are discussed here.

  • The different terms involved in the equations of motion of flexible bodies can be decomposed using shape integrals (see Appendix D3). Our framework readily supports this optional decomposition: it is the responsibility of the user to provide the terms and values of the expansion when numerical evaluation is to occur.

  • The definition of geometric stiffening requires attention in the general case. It is accounted for by the term kσ, presented in Appendix A. We discuss geometric stiffening in more detail in Appendix C.

  • The treatment of external loads was not addressed in detail in this article because the loads are application-specific (aerodynamics, hydrodynamics, etc.). The framework can accept external loads as arbitrary functions of multiple variables or as analytical expressions. In the former case, the user will have to provide an implementation of the function during the execution.

  • Even though the equations of motion are void of constraint forces, the values of these forces can be recovered. They can be expressed as functions of the external forces and the states of the system. It is not necessary to compute them by iteratively solving constraint equations.

  • The framework can easily include rheonomous constraints – for instance, for the pitch angle – without having to supply a dedicated torque. Pitch speed and accelerations can be directly introduced into the mechanical system if they are provided by a generic second-order pitch actuator model.

5.3 Limitations

In spite of the advantages listed in Sect. 5.1, the symbolic procedure presented in this work has some potential limitations.

  • Constraints and closed loops have currently not been added to the framework. The SymPy mechanics package supports additional constraint equations within Kane's method. We therefore hope that this limitation can be lifted in the future.

  • Large problems may challenge a symbolic calculation package: memory impact, calculation time, simplification times, and size of expressions may become significant. Some of these issues may be alleviated by introducing intermediate variables that are only substituted for in the numerical implementation or by using a recursive formulation of the solution procedure (Branlard2019).

We further note that the shape function approach is an approximate method: it introduces a separation of space and time early on in the development of the nonlinear equations of motion and applies low-order polynomial (usually linear or quadratic) approximations to eliminate high-order terms (see, e.g., Table 1 of Wallrapp1994). This was presented as an advantage in Sect. 5.1 because the equations are obtained in compact form and are readily linearized. Yet, the approximations introduced by the method may imply that nonlinearities are not well captured, which is why the models are labeled as “mid-fidelity” throughout this article. The domain of validity of the nonlinear or linear models presented may therefore be limited in time and space as opposed to fully nonlinear methods. Advanced methods to obtain high-fidelity reduced-order models from nonlinear dynamic systems are beyond the scope of this work; see, e.g., Steindl and Troger (2001), Benner et al. (2015), and Touzé et al. (2021).

6 Conclusions

We presented a symbolic framework to obtain the linear and nonlinear equations of motion of a multibody system made of rigid bodies, flexible bodies, and kinematic joints. Our approach is based on Kane's method and a nonlinear shape function representation of flexible bodies. We provided different wind energy examples and verified the results against OpenFAST simulations. The framework can readily provide models suitable to a wide range of applications with competitive computational times. The framework is open source, and the examples presented are available in the repository. Future work will focus on applying the framework to dedicated research projects, with more complex systems, and potentially extend the framework to account for closed-loop systems and arbitrary constraints.

Appendix A: Equations for a flexible body and shape integrals

In this section, we detail the equations of motion of a flexible body. The reader is referred to the following references for a complete treatment of the equations of motion: Shabana (2013), Schwertassek and Wallrapp (1999), and Wallrapp (1994). The subscript i, indicating the body index, is dropped. All quantities (vectors and matrices) are expressed in the body frame of reference; therefore, the prime notation is also dropped in this section. The number of flexible shape functions associated with the body is ne, the flexible degrees of freedom are qe, and the shape functions are gathered into a matrix Φ of size (3×ne). The equations of motion, given in Eq. (12), are repeated below:

(A1) M xx M x θ M xe M θ θ M θ e sym. M ee a i ω ˙ i q ¨ e + k ω , x k ω , θ k ω , e + 0 0 k e = f x f θ f e .

The different terms of the mass matrix are obtained as follows:


The integrals are volume integrals over the volume of the body (for beams, they reduce to line integrals). The notation [̃] represents the skew symmetric matrix. M is the mass of the body. The vector sCM is the vector from the origin of the body to undeflected center or mass (CM) of the body. The notations Ct (ne×3) and Cr (ne×3) are introduced to match Wallrapp's notations. The vector sP is the vector from the origin of the body to a deflected point of the body of elementary mass dm. The undeflected position of this point is written as sP0 and the displacement field u, such that sP=sP0+u. Typically, the displacement field is given by u=Φqe, but a higher-order expansion can also be introduced (see Wallrapp1994, and Appendix D4). Wallrapp also includes the elementary mass moment of inertia, which results in additional terms in the integrals (see Wallrapp1994). Such contributions are relevant, for instance, when considering the torsion of a beam (see Branlard2019). The block matrices Mxx, Mxe, and Mee do not depend on the deformation of the body and are therefore constant. The other terms are functions of qe. They may be expressed as linear combinations of constant integrals (see Appendix D3).

The quadratic velocity terms, kω, are given as




The first term of Eq. (A10) is obtained by vertically stacking the contribution of each shape function. In the standard input data format, this term is reshaped as the product OeΩ, where


The body elastic forces are given by

(A16) k e = k σ + K e q e + D e q ˙ e ,

where Ke and De are the elastic stiffness and damping matrices, and kσ represents geometric stiffening terms (see Appendix C). The elastic damping forces are often given as stiffness proportional damping. For more details, see Wallrapp (1994), and for more examples with elastic beams, see Branlard (2019). The external loads can be assumed to consist of distributed volume forces, p (in practice they are primarily surface forces or line forces), and a gravitational acceleration field, g. The components of the external loads in Eq. (A1) are then obtained by integration over the whole body:

Appendix B: Application of the shape function approach to an isolated beam

In this section, we illustrate how the elastic equations of Appendix A can be applied to an isolated beam. Examples of applications are further given in Sects. 4.3 and 4.4. We consider a beam directed along the z axis and bending in the x and y directions. Expressions are written in the coordinate system of the beam, and primes are dropped in this section. The beam properties are the following: length, L; mass per length, m; and bending stiffness, EIx and EIy. We assume that the displacement field is such that the shape functions are functions of z only: u(z,t)=i=1neΦi(z)qe,i(t). We also assume that the shape functions satisfy at least the geometric boundary conditions. The kinetic energy of the beam is T=120Lmu˙2dz=12ijMe,ijq˙e,jq˙e,i, where Me,ij is (see Eq. A7)

(B1) M e , i j = 0 L m ( z ) Φ i ( z ) Φ j ( z ) d z , i , j = 1 , n e .

Equation (B1) involves a scalar product of the shape functions at each spanwise position. Integrals over the moment of inertia can be used to account for torsion (see Branlard2019). The potential energy (strain energy) of the beam is obtained as V=12ijKe,ijqe,iqe,j, where Ke,ij represents the elements of the stiffness matrix, which, under the assumption of small deformations, are given by

(B2) K e , i j = 0 L [ EI y d 2 Φ i , x d z 2 d 2 Φ j , x d z 2 + EI x d 2 Φ i , y d z 2 d 2 Φ j , y d z 2 ] d z , i , j = 1 , n e .

Elongation and torsional strains (EA and GKt) can similarly be added to the strain energy and the stiffness matrix if longitudinal and torsional displacement fields are included in the shape functions. The external loads on the beam are assumed to consist of a distributed force vector, p(z). The virtual work done by the force p for each virtual displacement δqe,i provides the generalized force as (see Eq. A17)

(B3) f e , i = 0 L Φ i p d z .

The equations of motion of the isolated beam are then written in matrix form as

(B4) M e q ¨ e + D e q ˙ e + K e q e = f e ,

where qe=[qe,1,,qe,n]. Damping is typically added a posteriori to the equations, where the Rayleigh damping assumption is often used: De=αMe+βKe (stiffness proportional damping implies α=0). If the shape functions are mode shapes, then the shape functions are orthogonal, the mass and stiffness matrices are diagonal, and the stiffness values would be Ke,ii=ωe,i2Me,ii, with ωe,i=Ke,ii/Me,ii the eigenfrequency of the beam mode i. The modal damping is then given by De,ii=2ζiMe,iiωe,i, where ζi is the damping ratio associated with mode i.

If the beam is loaded axially by a force N(z) (assumed to be independent of the elastic degrees of freedom), then this force produces a distributed load in the transverse direction equal to n=zN(z)uz, with components in the y and z directions (see Branlard2019). The generalized force associated with this loading is then QN,i=0LΦindz. Inserting the expression of n and u, the generalized force has the form of a stiffness term: QN,i=-jKN,ijqe,j, with

(B5) K N , i j = - 0 L Φ i d d z N ( z ) d Φ j d z d z = 0 L N ( z ) d Φ i d z d Φ j d z - N ( z ) Φ i d Φ j d z 0 L ,

where integration by parts was used to obtain the second equality. Examples of applications are given in Sects. 4.3 and 4.4. The fact that an axial load leads to a stiffness term is referred to as “geometric stiffness”, which is the topic of Appendix C.

Appendix C: Geometric stiffness

C1 General treatment

Geometric stiffness refers to the apparent change of stiffness of a structure depending on the loading it is subject to. In this section, we present a linear formulation of geometric stiffness for a flexible body undergoing motion and subject to arbitrary loading, inspired by Schwertassek and Wallrapp (1999). Additional details may be found in Wallrapp and Schwertassek (1991). The main component of the geometric stiffening term kσ can be written as

(C1) k σ = K g q e ,

where Kg is the geometric stiffness matrix of shape ne×ne. In general, this matrix is time-dependent, as it is a function of the inertial and external loads acting on the body. The inertial loads consist of contributions from the linear acceleration (a), rotational acceleration (ω˙), and cross products of the rotational velocity of the body (centrifugal and gyroscopic terms). The external loads consist of the gravitational force, distributed forces per unit length (p), point loads (Fk), and point moments (τk), where k is the node index where the point loads are applied. Each of these contributions can be computed at each time step using a linear superposition of unit geometric stiffness matrices, denoted Kg, as follows:

(C2) K g = α = 1 3 a α - g α K gt , α + ω ˙ α K gr , α + α = 1 3 β = 1 3 ω α ω β K g ω , α β + α = 1 3 p α K gp , α + k F α k K gF , α k + τ α k K g τ , α k ,

where the indices α and β run on the x, y, and z coordinates of the body reference frame. The matrices Kg,α or Kg,αβ have the shape ne×ne and are obtained as the geometric stiffness matrices for unit accelerations, loads, or products of rotational velocities in the given direction defined by α and β (x, y, or z). For instance, Kgt,z is the geometric stiffness matrix corresponding to a unit acceleration in the z direction, Kgω,xyk is the geometric stiffness matrix corresponding to a unit gyration about the x and y directions (centrifugal effect), and KgF,xk is the geometric stiffness matrix corresponding to a unit force in the x direction applied at the node k along the body. The effect of the Coriolis force is not mentioned in the work of Schwertassek and Wallrapp and not explicitly accounted for in Eq. (C2). The Coriolis force, 2m(z)ω×(jΦj(z)q˙e,j), is proportional to q˙e. Because the instantaneous beam slope is proportional to qe, the geometric stiffening term consists of terms of the form qe,jq˙e,k. In Table 1 of Wallrapp (1994), it is stated that nonlinear terms of the form qe,jq˙e,k are neglected. Yet, if the steady-state deflection qe is significant, then the influence of the Coriolis term on the geometric stiffening may be significant. The effect can be included as an additional term in Eq. (C1) that is a function of qe and q˙e (expressions are provided in Eq. C2). We note that the terms Kg have different units; for instance, the terms Kgt, are expressed in N s2 m−2 .

C2 Expressions for a beam directed along z

The expression for each of these matrices are given in Schwertassek and Wallrapp (1999) in the context of the finite-element method. The general expressions for a shape function approach would be beyond the scope of this article, but we provide the expressions for a beam below.

We adopt the same notations as Appendix B to describe the flexible beam. Following the developments that led to Eq. (B5), the geometric correction associated with an axial load N is given by the generalized force:

(C3) k σ , N , i = 0 L N ( z ) Φ i j d Φ j d z q j d z .

When the axial load is not a function of the degrees of freedom, this expression can be expressed as a stiffness matrix, as indicated in Eq. (B5). The different unit geometric matrices introduced in Appendix C can be determined using a form of Eq. (B5), where the axial load N is replaced by the unit inertial or external load. Since the beam is directed along the z direction, we focus on the terms where the loads act in the z direction, with all other terms being zero or negligible. The ij component of the matrix Kgt,z is obtained by considering a unit vertical acceleration:

(C4) K gt , z , i j = 0 L N ( z ) d Φ i d z d Φ j d z d z , N ( z ) = z L m ( z ) d z .

We write zk as the coordinate of node k along the beam. The ij component of the matrix KgF,zk is obtained as

(C5) K gF , z , i j k = 0 L N ( z ) d Φ i d z d Φ j d z d z , N ( z ) = 1  if  z < z k , 0  otherwise .

The ij component of the matrix Kgω,αβ is obtained by considering unit centrifugal loads generated using independent rotations around the unit vectors ex, ey, and ez:

(C6) K g ω , α β , i j = 0 L - e z e ̃ α e ̃ β N ( z ) d Φ i d z d Φ j d z d z , N ( z ) = z L m ( z ) s P 0 d z .

Similarly, the ij component of the matrix Kgr,α is

(C7) K gr , α , i j = 0 L - e z e ̃ α N ( z ) d Φ i d z d Φ j d z d z , N ( z ) = z L m ( z ) s P 0 d z .

The Coriolis force, 2m(z)ω×(jΦj(z)q˙e,j), can also have an axial contribution. Using Eq. (C3), the generalized force is

(C8) k σ , Cor , i = 0 L N ( z ) Φ i j d Φ j d z q j d z , N ( z ) = 2 z L m ( z ) k ω x Φ k , y ( z ) - ω y Φ k , x ( z ) q ˙ k d z .

Equation (C8) may be rearranged by introducing a three-dimensional tensor made of shape integrals that are independent of the elastic degrees of freedom and the rotational speed and therefore speed up the evaluation of this expression at each time step. The vector kσ,Cor(qe,q˙e) is added to the right-hand side of Eq. (C1).

C3 Integration into the equations of motion

The term kσ=Kgqe appears on the third block row of the equations of motion of the flexible body (Eq. A1). Because of the linearity with respect to the acceleration, rotational velocities, and forces, the different contributions can optionally be incorporated into the third block row of the mass matrix (Me), the term kω,e, and the term fe, respectively. For instance, the term aαKgt,αqe can be reorganized as [Kgt]qea (using loose notations); therefore, the mass matrix can be updated such that Mxe becomes Mxe+[Kgt]qe. When a Taylor expansion is used, such integration is easily implemented as a first-order term (see Appendix D3).

Appendix D: Alternative formulations

Different formulations of flexible multibody dynamics using shape functions are found in the literature. Some of the alternatives are briefly discussed in this section.

D1 Jacobian and velocity transformation matrix

In Eq. (7), the Jacobian terms J and the virtual work are expressed in vector form. In such form, there is no need to state in which coordinate system the different vectors are expressed. This is convenient to reduce the size of the expressions when using symbolic calculations. In a numerical framework, the vector will have to be expressed in a common frame. When such an approach is used (see, e.g., Lemmer2018, and Branlard2019), the Jacobians are sometimes stacked into a matrix form:

(D1) J = J v J ω J e .

Some implementation choices are needed depending on whether these matrices are expressed in the global frame or a body frame. The Jacobian matrices are referred to as “velocity transformation matrix”, and the link between formulations in global and local coordinates is given in Branlard (2019). In the same reference, recursive relationships are given for tree-like assembly of bodies to help express the Jacobian matrices of each body recursively, based on the matrices of the parent body. It is also noted that the quadratic velocity terms, kω, can be obtained using the time derivative of the Jacobian matrix.

D2 Rotations and torsion

In this article, we have not elaborated on the change of orientation introduced by shape functions. In most applications, bodies are connected at their extremities, and the deflection slope at a body extremity will induce a rotation of the subsequent body (e.g., tilting and rolling of the nacelle at the tower top). The deflection slope can be obtained form the knowledge of the shape functions. This is readily accounted for by introducing a time-varying rotation matrix between bodies, and this is the approach used in our symbolic framework. A formalism of rotations of bodies connected at their extremities is given in Branlard (2019). A more general formulation, introducing shape function rotations Ψ, is given in (Wallrapp1994; Schwertassek and Wallrapp1999; Lemmer2018). In such a formulation, the linear rotation field is obtained as I+Ψq̃, where I is the identity matrix.

D3 Shape integrals and Taylor expansion

The results presented in Appendix A consist of integrals over the displaced points of the structure, sP=sP0+u, where the displacement field is u=Φqe. The undeflected position of the structure (sP0) is constant, and the shape functions are known at the initialization; the only time-varying terms are the degrees of freedom qe. Therefore, the integrals can be precomputed by decomposing them into a constant part and a part that is linear with respect to the degrees of freedom qe. The precomputed integrals are referred to as “shape integrals”. For a given term T (standing, for instance, for Mθ,θ, Ct, Cr, Gr, Ge, or Oe), the shape integral expansion is

(D2) T ( q e ) = T 0 + j = 1 . . n e T j 1 q e , j .

If T is an array, T0 and Tj1 have the same shape as T. As an example, the application of the shape integral expansion to the term Mxθ (see Eq. A3) gives

(D3) M x θ = - s ̃ P d m = M x θ 0 + j = 1 . . n e M x θ , j 1 q e , j ,


(D4) M x θ 0 = - s ̃ P 0 d m , M x θ , j 1 = - Φ ̃ j d m .

The zeroth- and first-order shape integrals always consist of integrals over the components of sP0 and Φ, which can be precomputed for a given flexible body. We note that the precomputed shape integrals can in turn be obtained from intermediate integrals (e.g., the S and N terms introduced by Wallrapp1994, or the σ, Σ, Υ, and Ψ terms introduced by Shabana2013). The zeroth- and first-order shape integrals are stored using a Taylor object-oriented class in the standard input data format defined by Wallrapp. The YAMS library can compute the shape integrals using a direct integration or using a finite-element formulation (see Schwertassek and Wallrapp1999).

The geometric stiffness introduced in Appendix C is linear in the elastic degrees of freedom qe. Therefore, the unit geometric stiffness matrices (which are also shape integrals) can be conveniently added into the first-order terms of Eq. (D2). For instance, if we write Mex (given in Eq. A6) using a first-order expansion, Mex=Mex0+Mex1qe, then the geometric stiffening effect can directly be inserted into the first-order term, such that Mex1 becomes Mex1+Kgt. Similarly, the term Kgr can be inserted into Mθe1, Kgω into Oe1, KgF into Φ1, and Kgτ into Ψ1 in the calculation of the generalized forces. The different contributions are summarized in Table 6.9 of the book of Schwertassek and Wallrapp (1999). A shortcoming of inserting the geometric stiffness effects into the first-order coefficient is that it could make the mass matrix symmetric (if the user code assumes Mxe=Mext), instead of acting only on the third block row of the mass matrix.

D4 Taylor expansion of the displacement field

In the work of Wallrapp (Wallrapp1993, 1994), the displacement field is assumed to be a function of the degrees of freedom, u=Φu(qe)qe, where Φu consists of a Taylor series expansion of the shape functions that contain Φ0 and Φ1 terms. The resulting equations of motion are still expressed using shape integrals of the form given in Eq. (D2), but the 1 terms will contain some additional integrals over Φ1. The advantage of this method is that the Φ1 terms effectively account for the geometric stiffness. In practice, it is equivalent, and as convenient, to neglect the Φ1 terms and introduce the geometric stiffness using the method presented in Appendix C (and optionally integrate them into the 1 terms as presented in Appendix D3).

D5 ElastoDyn and the partial loads approach

The ElastoDyn module of OpenFAST (Jonkman et al.2021) uses the so-called “partial loads” approach to implement the equations of motion. The underlying theory used to derive the equations of motion is the same as Kane's formalism presented in Sect. 2, but the partial load approach takes advantage of the fact that the calculation of reaction loads or point loads at body extremities requires similar terms to the ones needed for the equations of motion. In the discussion below, we assume that the different bodies of the structure form a tree structure with the root at the bottom and the leaves above. For a tree-like structure, there is a natural relationship between loads in the structure and the degrees of freedom. A virtual displacement of a given degree of freedom will only displace the structure above it. The equation of motion of this degree of freedom can therefore be obtained from the virtual work of the loads at a point located just above the degree of freedom, as if the entire structure above was replaced by lumped loads. The point loads contain contributions from the external loads above the point in consideration but also inertial and gyroscopic loads associated with all the degrees of freedom of the system. If the point is at a joint, the loads corresponds to the reaction loads at this point. We write P as the point located after a given degree of freedom r. The equation of motion for this degree of freedom is obtained as if the system was isolated:

(D5) f r + f r = 0 = J v P , r f P + J ω P , r τ P + h r ,

where JvP,r and JωP,r are the partial velocities of point P with respect to the degree of freedom r, fP and τP are three vectors containing the force and torque from the structure above the degree of freedom r (including external and inertial contributions), and hr is the generalized load associated with the isolated degree of freedom r (e.g., the elastic loads for a flexible body, or the spring and damping loads for a degree of freedom representing a joint). The point loads fP and τP can be decomposed into terms that are proportional to the accelerations of all the degrees of freedom (indexed with r) and additional terms (labeled “t”):

(D6) f P = j = 1 n q f P , j q ¨ j + f P , t , τ P = j = 1 n q τ P , j q ¨ j + τ P , t .

The terms fP,r and τP,r act as generalized masses, and they are referred to as “partial loads”. Combining Eqs. (D5) and (D6), the term rj of the mass matrix and the term r of the right-hand side of the equation of motion (Eq. 23) are obtained as

(D7) M rj = - J v P , r f P , j - J ω P , r τ P , j , F r = J v P , r f P , t + J ω P , r τ P , t + h r .

Therefore, the knowledge of the partial loads and the partial velocities at key points of the structure (typically, points where user outputs are desired) can be used to obtain the reaction loads (Eq. D6) and the equations of motion (Eq. D7). This is the approach used in ElastoDyn: the loads at key points of the structure were derived using hand calculations, and then the partial loads were used for the implementation of the outputs and the equations of motion. The reader is referred to the notes provided in the online documentation of ElastoDyn for more details (Jonkman et al.2021). A general procedure to obtain partial loads can be devised (using kinematics to find velocities and acceleration in the structure and computing the loads from the tree top to the root) but would be beyond the scope of this article.

Appendix E: Equations of motion of simple wind turbine models

In this section, we present the equations of motion for the examples presented in Sect. 4.

E1 Two-degrees-of-freedom model of a land-based or fixed-bottom wind turbine

In this section, we provide some intermediate values to obtain the equations of motion given in Sect. 4.4. We use the hat notation to indicate unit vectors of a frame, where the frame is identified as t, n, and r for the tower, nacelle, and rotor, respectively. For instance, vt^x is the unit vector in the x direction of the tower frame. The degrees of freedom are q=(q,ψ). The kinematics of the tower (at its origin) are zero:

(E1) v O , T = 0 , ω T = 0 , a O , T = 0 .

All Jacobians are zero except Je,1T=1. The inertial force, torque, and elastic force are

(E2) f T = C t T x q ¨ t ^ x + M T g t ^ z , τ T = C r T y q ¨ t ^ y , E T = f e + D e q ˙ + ( K e + K q ) q + M e q ¨ .

The nacelle kinematics (at its center of mass) are


The Jacobians with respect to q are

(E5) J v , 1 N = t ^ x + ν y z NG n ^ x - ν y x NG n ^ z , J ω , 1 N = ν y t ^ y .

The inertial force and torque on the nacelle are

(E6) f N = M N q ¨ t ^ x + M N - ν y 2 x NG q ˙ 2 + ν y z NG q ¨ n ^ x + M N - ν y 2 z NG q ˙ 2 - ν y x NG q ¨ n ^ z , τ N = J y,N ν y q ¨ n ^ y .

The kinematics of the rotor are


The corresponding Jacobians with respect to q (“1”) and ψ (“2”) are


The inertial force and torque on the rotor are


E2 Three-degrees-of-freedom model of a land-based or fixed-bottom wind turbine

The equations of motion for the model presented in Sect. 4.5, with q=(q1,q2,ψ), are given in this section. The elements of the mass matrix are


The elements of the forcing vector are


E3 Three-degrees-of-freedom model of a floating wind turbine

The equations of motion for the model presented in Sect. 4.6, with q=(x,ϕ,qT), are given in this section. The elements of the mass matrix are


The elements of the forcing vector are

Code availability

The latest source code of YAMS is available on GitHub as a subpackage of the Wind Energy Library, WELIB (, Branlard2022a). A static version is available on Zenodo (, Branlard2022b). The examples given in this articles are found in the folder welib/yams/papers of the repository.

Author contributions

Both authors exchanged ideas over the last 2 years about the implementation of such a framework and its application to wind energy. EB wrote a Python implementation and JG wrote a Maxima implementation. EB wrote the main part of the article, with feedback and contributions from JG.

Competing interests

The contact author has declared that none of the authors has any competing interests.


The views expressed in the article do not necessarily represent the views of the DOE or the U.S. Government.

Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.


Funding was provided by the U.S. Department of Energy Office of Energy Efficiency and Renewable Energy Wind Energy Technologies Office.

Financial support

This work was funded under the Technology Commercialization Fund Project, supported by the DOE's Wind Energy Technologies Office.

Review statement

This paper was edited by Raimund Rolfes and reviewed by two anonymous referees.


ANSYS: (last access: 19 March 2022), 2022. a

Bauchau, O. A.: Flexible Multibody Dynamics, Solid Mechanics and Its Applications, Springer, Dordrecht,, 2011. a

Benner, P., Gugercin, S., and Willcox, K.: A Survey of Projection-Based Model Reduction Methods for Parametric Dynamical Systems, SIAM Rev., 57, 483–531,, 2015. a

Bielawa, R.: Rotary wing structural dynamics and aeroelasticity, AIAA education series, American Institute of Aeronautics and Astronautics, 1801 Alexander Bell Drive, Reston, VA 20191, 2006. a

Branlard, E.: Flexible multibody dynamics using joint coordinates and the Rayleigh-Ritz approximation: The general framework behind and beyond Flex, Wind Energy, 22, 877–893,, 2019. a, b, c, d, e, f, g, h, i, j, k, l, m

Branlard, E.: WELIB, Wind Energy Library, GitHub repository (last access: November 2022), 2022a. a, b, c

Branlard, E.: WELIB, Zenodo,, 2022b. a

Branlard, E., Giardina, D., and Brown, C. S. D.: Augmented Kalman filter with a reduced mechanical model to estimate tower loads on a land-based wind turbine: a step towards digital-twin simulations, Wind Energ. Sci., 5, 1155–1167,, 2020a. a

Branlard, E., Jonkman, J., Dana, S., and Doubrawa, P.: A digital twin based on OpenFAST linearizations for real-time load and fatigue estimation of land-based turbines, J. Phys. Conf. Ser., 1618, 022030,, 2020b. a

Docquier, N., Poncelet, A., and Fisette, P.: ROBOTRAN: a powerful symbolic gnerator of multibody models, Mech. Sci., 4, 199–219,, 2013. a

Gede, G., Peterson, D., Nanjangud, A., Moore, J., and Hubbard, M.: Constrained Multibody Dynamics With Python: From Symbolic Equation Generation to Publication, in: Proceedings of the ASME 2013 International Design Engineering Technical Conferences and Computers and Information in Engineering Conference, Portland, Oregon, USA, 4–7 August 2013,, 2013. a, b, c

Geisler, J.: CADynTub: Wind Turbine Model from OpenFAST Data using CADyn Equations of Motion, Github, (last access: November 2022), 2021. a

Géradin, M. and Cardona, A.: Flexible Multibody Dynamics: A Finite Element Approach, Wiley, New York, ISBN 978-0-471-48990-0, 2001. a

Jelenić, G. and Crisfield, M.: Geometrically exact 3D beam theory: implementation of a strain-invariant finite element for statics and dynamics, Comput. Method. Appl. M., 171, 141–171,, 1999. a

Jonkman, B., Mudafort, R. M., Platt, A., Branlard, E., Sprague, M., Jonkman, J., Vijayakumar, G., Buhl, M., Ross, H., Bortolotti, P., Masciola, M., Ananthan, S., Schmidt, M. J., Rood, J., Damiani, R., Mendoza, N., Hall, M., and Corniglion, R.: OpenFAST v3.1.0. Open-source wind turbine simulation tool, Zenodo,, 2021. a, b, c, d, e

Jonkman, J., Butterfield, S., Musial, W., and Scott, G.: Definition of a 5 MW Reference Wind Turbine for Offshore System Development, Tech. Rep. NREL/TP-500-38060, National Renewable Energy Laboratory,, 2009. a, b

Jonkman, J. M.: Dynamics of offshore floating wind turbines–model development and verification, Wind Energy, 12, 459–492,, 2009. a

Jonkman, J. M., Branlard, E. S. P., and Jasa, J. P.: Influence of wind turbine design parameters on linearized physics-based models in OpenFAST, Wind Energ. Sci., 7, 559–571,, 2022. a, b

Kane, T. R. and Wang, C. F.: On the Derivation of Equations of Motion, J. Soc. Ind. Appl. Math., 13, 487–492,, 1965. a, b

Kurz, T. and Eberhard, P.: Symbolic Modeling and Analysis of Elastic Multibody Systems, in: International Symposium on Coupled Methods in Numerical Dynamics Split, Croatia, September 16-19, 2009. a, b

Kurtz, T., Eberhard, P., Henninger, C., and Schiehlen, W.: From Neweul to Neweul-M2: symbolical equations of motion for multibody system analysis and synthesis, Multibody Syst.Dyn., 24, 25–41,, 2010. a

Lange, C., Kövecses, J., and Gonthier, Y.: Benchmarking of Multibody System Simulations: Points to Consider, in: CcToMM Symposium on Mechanisms, Machines, and Mechatronics, Saint-Hubert, Quebéc, 31 May–1 June 2007. a

Lemmer, F.: Low-order modeling, controller design and optimization of floating offshore wind turbines, PhD thesis, Universität Stuttgart, (last access: November 2022), 2018. a, b, c, d

MBDyn: (last access: 19 March 2022), 2022. a

Merz, K. O.: STAS Aeroelastic 1.0 – Theory Manual., Tech. rep., SINTEF Energi AS., Trondheim, 2018. a

MotionGenesis: MotionGenesisTM Kane Tutorial, Tech. rep., Motion Genesis LLC, (last access: November 2022), 2016. a

Øye, S.: Fix Dynamisk, aeroelastisk beregning af vindmøllevinger, Report AFM83-08, Fluid Mechanics, DTU, Lyngby, Denmark, 1983. a

Pöschke, F., Gauterin, E., Kühn, M., Fortmann, J., and Schulte, H.: Load mitigation and power tracking capability for wind turbines using linear matrix inequality-based control design, Wind Energy, 23, 1792–1809,, 2020. a

Reckdahl, K. and Mitiguy, P.: Autolev Tutorial, Tech. rep., OnLine Dynamics Inc., Sunnyvale, CA, 1996. a

Schwertassek, R. and Wallrapp, O.: Dynamik flexibler Mehrkörpersysteme, Friedr. Vieweg & Sohn, Braunschweig, 1999 (in German). a, b, c, d, e, f, g, h

Shabana, A.: Dynamics of Multibody Systems, Dynamics of Multibody Systems, Cambridge University Press, Cambridge CB2 8BS, United Kingdom, ISBN 9781107042650, 2013. a, b, c, d, e

Simani, S.: Advanced Issues of Wind Turbine Modelling and Control, J. Phys. Conf. Ser., 659, 012001,, 2015. a

Simo, J.: A finite strain beam formulation. The three-dimensional dynamic problem. Part I, Comput. Method. Appl. M., 49, 55–70,, 1985. a

SIMPACK: (last access: 19 March 2022), 2022. a

Sønderby, I. and Hansen, M. H.: Open-loop frequency response analysis of a wind turbine using a high-order linear aeroelastic model, Wind Energy, 17, 1147–1167,, 2014. a

Steindl, A. and Troger, H.: Methods for dimension reduction and their application in nonlinear dynamics, Int. J. Solids Struct., 38, 2131–2147,, 2001. a

SymPy: (last access: November 2022), 2021.  a, b

Touzé, C., Vizzaccaro, A., and Thomas, O.: Model order reduction methods for geometrically nonlinear structures: a review of nonlinear techniques, Nonlinear Dynam., 105, 1141–1190,, 2021. a

Verlinden, O., Kouroussis, G., and Conti, C.: EasyDyn: a framework based on free symbolic and numerical tools for teaching multibody systems, in: Multibody Dynamics 2005, ECCOMAS Thematic Conference, Madrid, Spain, 21–24 June 2005. a

Wallrapp, O.: Standard Input Data of Flexible Members in Multibody Systems, in: Advanced Multibody System Dynamics. Solid Mechanics and Its Applications, edited by: Schiehlen, W., vol. 20, pp. 445–450, Springer, Dordrecht,, 1993. a

Wallrapp, O.: Standardization of flexible body modeling in multibody system codes, part i: Definition of standard input data, Journal of Structural Mechanics, 22, 283–304, 1994. a, b, c, d, e, f, g, h, i, j, k, l

Wallrapp, O. and Schwertassek, R.: Representation of geometric stiffening in multibody system simulation, Int. J. Numer. Meth. Eng., 32, 1833–1850,, 10.1002/(ISSN)1097-0207, 1991. a


We address the first-order approximation in Appendix D4.


The relevant equations of the shape function approach for a beam are given in Appendix B.

Short summary
The article presents a framework to obtain the linear and nonlinear equations of motion of a multibody system including rigid and flexible bodies. The method yields compact symbolic equations of motion. The applications are many, such as time-domain simulation, stability analyses, frequency domain analyses, advanced controller design, state observers, and digital twins.
Final-revised paper