A symbolic framework to obtain midfidelity models of flexible multibody systems with application to horizontalaxis wind turbines
 ^{1}National Renewable Energy Laboratory, Golden, CO 80401, USA
 ^{2}Energy and Biotechnology, Flensburg University of Applied Sciences, 24943 Flensburg, Germany
 ^{1}National Renewable Energy Laboratory, Golden, CO 80401, USA
 ^{2}Energy and Biotechnology, Flensburg University of Applied Sciences, 24943 Flensburg, Germany
Correspondence: Emmanuel Branlard (emmanuel.branlard@nrel.gov)
Hide author detailsCorrespondence: Emmanuel Branlard (emmanuel.branlard@nrel.gov)
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 midfidelity 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 tradeoff to obtain midfidelity 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 landbased 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 opensource project.
This work was authored in part by the National Renewable Energy Laboratory, operated by Alliance for Sustainable Energy, LLC, for the U.S. Department of Energy (DOE) under contract no. DEAC3608GO28308. The U.S. Government retains and the publisher, by accepting the article for publication, acknowledges that the U.S. Government retains a nonexclusive, paidup, irrevocable, worldwide license to publish or reproduce the published form of this work, or allow others to do so, for U.S. Government purposes.
The next generation of wind turbine digital technologies and control systems require versatile aeroservohydroelastic 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, highperformance or embedded code (for standalone 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 (Øye, 1983; Branlard, 2019) 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 reducedorder models from these tools or extend the models to additional degrees of freedom.
Tools with linearization capabilities, such as HAWCStab2 (Sønderby and Hansen, 2014) or OpenFAST (Jonkman et al., 2021), are dedicated to horizontalaxis wind turbines, and the evaluation of the gradients are limited to hardcoded analytical expressions or numerical finite differences. Small implementation changes often require extensive redevelopment, and the range of applications of the tools remains limited (Simani, 2015). 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 codesign 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 Wang, 1965) and a nonlinear shape function representation of flexible bodies (Shabana, 2013) described using a standard input data (SID) format (Wallrapp, 1994; Schwertassek and Wallrapp, 1999). 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 codesign approaches. We implemented an opensource version in Python using SymPy (SymPy, 2021), leveraging its mechanical toolbox. Alternative symbolic frameworks found in the literature are usually limited to rigid bodies (Verlinden et al., 2005; Kurz and Eberhard, 2009; Gede et al., 2013; Docquier et al., 2013) or are closedsource or using proprietary software (Reckdahl and Mitiguy, 1996; Kurtz et al., 2010; MotionGenesis, 2016; Lemmer, 2018).
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 (Simo, 1985; Jelenić and Crisfield, 1999; Géradin and Cardona, 2001; Bauchau, 2011) 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 (ANSYS, 2022), SIMPACK (SIMPACK, 2022), or MBDyn (MBDyn, 2022). These more advanced approaches target different applications than those envisioned in this study: they are suitable for numerical simulations, but they cannot provide symbolic midfidelity 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.
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 n_{b} 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 n_{q}, 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., Branlard, 2019). 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 $\mathit{q},\dot{\mathit{q}}$, and $\ddot{\mathit{q}}$, where $\left(\phantom{\rule{0.125em}{0ex}}\dot{\phantom{\rule{0.125em}{0ex}}}\phantom{\rule{0.125em}{0ex}}\right)$ 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., Shabana, 2013)
where r_{i}, v_{i}, and a_{i} are the position, velocity, and acceleration of the origin of the body, respectively; ${\mathit{s}}_{{P}_{\mathrm{0}}}$ 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 P_{0} for the undeformed position; u_{P} 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 $\left(\dot{\phantom{\rule{0.33em}{0ex}}}\right)$ and $(\dot{\phantom{\rule{0.33em}{0ex}}}{)}_{\mathrm{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 R_{i} that transforms coordinates from the body frame to the inertial frame, and by definition $\left[{\stackrel{\mathrm{\u0303}}{\mathit{\omega}}}_{\mathrm{i}}\right]={\dot{\mathit{R}}}_{\mathrm{i}}{\mathit{R}}_{\mathrm{i}}^{\mathrm{T}}$, where $\left[\stackrel{\mathrm{\u0303}}{\phantom{\rule{0.33em}{0ex}}\phantom{\rule{0.33em}{0ex}}}\right]$ 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 matrixvector 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 Wang, 1965) is a powerful and systematic way to obtain the equations of motion of a system. The procedure leads to n_{q} coupled equations of motion:
where ${f}_{r}^{\ast}$ is associated with inertial loads and f_{r} 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:
The terms f_{ri} and ${f}_{\text{ri}}^{\ast}$ 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 f_{ri} and ${f}_{\text{ri}}^{\ast}$. The inertial force, ${\mathit{f}}_{\mathrm{i}}^{\ast}$, and inertial torque, ${\mathit{\tau}}_{\mathrm{i}}^{\ast}$, acting on the body are
where m_{i} is the mass of the body, a_{G,i} is the acceleration of its center of mass with respect to the inertial frame, and I_{G,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 ${f}_{\text{ri}}^{\ast}$ is defined as
with
where v_{G,i} is the velocity of the body mass center with respect to the inertial frame. The partial velocities, or Jacobians, J_{v} 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 ${f}_{\text{ri}}^{\ast}$ 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 f_{i} and τ_{i}. The component f_{ri} is then given by
Equivalently, the contributions from each individual force, f_{i,j}, acting on a point P_{j} of the body i, and each torque, τ_{i,k}, can be summed using the appropriate partial velocity to obtain f_{ri}:
where ${\mathit{v}}_{{P}_{j}}$ 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 f_{ri} and ${f}_{\text{ri}}^{\ast}$. 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 O_{i}. 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 zerothorder^{1} Rayleigh–Ritz approximation, the displacement field at a given point, P, is given by the sum of shape function contributions:
where Φ_{ij} are the shape functions (displacement fields) of body i, and q_{e,ij} is the subset of q consisting of the elastic coordinates of body i, of size n_{e,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 (Wallrapp, 1994)
where the x, θ, and e subscripts respectively indicate the translation, rotation, and elastic components; M is the mass matrix of dimension $\mathrm{6}+{n}_{\mathrm{e},i}$ made of the block matrices ${\mathit{M}}_{\text{xx}},\mathrm{\dots},{\mathit{M}}_{\text{ee}}$; a_{i} and ${\dot{\mathit{\omega}}}_{\mathrm{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; k_{e} 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_{ω}, k_{e}, and f are given in Appendix A. These terms depend on q, $\dot{\mathit{q}}$, and Φ_{i}. The inertial force, torque, and elastic loads are
The external and elastic loads are
The components of ${f}_{\text{ri}}^{\ast}$ and f_{ri}, for r=1…n_{q}, are then defined as
with
where v_{O,i} is the velocity of the body with respect to the inertial frame. The term J_{e,ri} consists of 0 and 1 because q_{e,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 n_{q} equations of motion given in Eq. (4) are gathered into a vertical vector e. They are recast into the form
or
where $\mathit{M}=\frac{\partial \mathbf{e}}{\partial \ddot{\mathit{q}}}$ is the system mass matrix and F is the forcing term vector – that is, the remainder terms of the equation ($\mathit{F}=\mathbf{e}+\mathit{M}\ddot{\mathit{q}})$. In this section, we introduce the vector u to represent the timedependent 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 $\dot{\mathit{q}}$, respectively. The nonlinear equation given in Eq. (23) is easily integrated numerically, for instance by recasting the system into a firstorder system or by using a dedicated secondorder system time integrator.
In various applications, a linear timeinvariant 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
The linearized equations about this operating point are obtained using a Taylor series expansion:
with
where M_{0}, C_{0}, and K_{0} are the linear mass, damping, and stiffness matrices, respectively; Q_{0}δu is the linear forcing vector (Q_{0} 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 (${\ddot{\mathit{q}}}_{\mathrm{0}}=\mathrm{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 parametervarying models used in many advanced control applications.
In this section, we discuss the Python opensource symbolic calculation framework that we implemented according to the equations given in Sect. 2. A Maxima implementation from the same authors is also available (Geisler, 2021).
The Python library YAMS (Yet Another Multibody Solver) started as a numerical tool published in previous work (Branlard, 2019). 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 (SymPy, 2021). 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 firstorder 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.

The lowlevel 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.

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.

The third level consists of template models such as generic landbased 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.
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 Pythonready 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 $\dot{\mathit{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 finiteelement 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 timeintegrate 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 (Branlard, 2022a). The repository contains tests and working examples, including the ones presented in Sect. 4.
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 hydroforce 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 (Branlard, 2022a).
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; M_{T}, mass; L_{T}, length; $({J}_{x,\mathrm{T}},{J}_{y,\mathrm{T}},{J}_{z,\mathrm{T}})$, diagonal coefficients of the inertia tensor about the center of gravity and in body coordinates; r_{TG}, vector from body origin to body center of mass of coordinates $({x}_{\text{TG}},{y}_{\text{TG}},{z}_{\text{TG}})$ 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 L_{B}=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): ${\mathbf{\Phi}}_{\mathrm{1}}={\mathrm{\Phi}}_{\mathrm{1},\mathrm{x}}{\mathit{e}}_{\mathrm{x}}=\mathrm{\Phi}{\mathit{e}}_{\mathrm{x}}$, where e_{x} 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 $\mathrm{\Phi}\left(r\right)=(r/R{)}^{\mathrm{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): $\mathit{u}(r,t)={\mathbf{\Phi}}_{\mathrm{1}}\left(r\right){q}_{\mathrm{e},\mathrm{1}}\left(t\right)=\mathrm{\Phi}\left(r\right)q\left(t\right){\mathit{e}}_{\mathrm{x}}$, where q_{e,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 p_{x}(r). The generalized mass and stiffness are computed based on the mass per length (m) and flapwise bending stiffness (EI_{y}) of the blade, according to Eq. (B1):
The generalized force is obtained from Eq. (B3):
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:
The geometric stiffness contribution of the axial load is obtained from Eq. (B5) as
The geometric stiffness, K_{g}, 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:
where k_{Ω} is referred to as the “rise factor” or “Southwell coefficient”, and in our approximation, it is found to be constant: ${k}_{\mathrm{\Omega}}={K}_{\mathrm{g}}\left(\mathrm{\Omega}\right)/{M}_{\mathrm{e}}/{\mathrm{\Omega}}^{\mathrm{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., $(\mathrm{\Omega}/\mathit{\omega}{)}^{\mathrm{2}}\lesssim \mathrm{5}$; see Bielawa, 2006), which is the case for wind energy applications.
The treatment for a shape function purely in the blade edgewise direction is similar, using Φ_{2}=Φ_{2}e_{θ}, 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}_{\mathit{\theta},\text{centri}}\left(r\right)={\mathrm{\Omega}}^{\mathrm{2}}{u}_{\mathit{\theta}}\left(r\right)\mathrm{d}m\left(r\right)$, with u_{θ}=Φ_{2}q_{2}. This leads to a generalized force equal to ${\int}_{\mathrm{0}}^{L}{p}_{\mathit{\theta},\text{centri}}{\mathrm{\Phi}}_{\mathrm{2}}\mathrm{d}r={\mathrm{\Omega}}^{\mathrm{2}}{M}_{\mathrm{e}}{q}_{\mathrm{2}}$, or, equivalently, to a stiffness term: ${K}_{\mathit{\omega}}={\mathrm{\Omega}}^{\mathrm{2}}{M}_{\mathrm{e}}$. It can be verified that this generalized force corresponds to the contribution ${O}_{\mathrm{e},\mathrm{11}}{\mathit{\omega}}_{\mathrm{x}}^{\mathrm{2}}$, from k_{ω,e}, given in Eq. (A10). For an edgewise mode, the frequency therefore evolves as
with ${k}_{\mathrm{\Omega}}={K}_{\mathrm{g}}\left(\mathrm{\Omega}\right)/{M}_{\mathrm{e}}/{\mathrm{\Omega}}^{\mathrm{2}}$ and with K_{g} 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 $\text{diag}\left({\mathit{M}}_{\mathrm{e}}\right)=[\mathrm{9.5}\times {\mathrm{10}}^{\mathrm{3}},\mathrm{1.5}\times {\mathrm{10}}^{\mathrm{4}},\mathrm{5.7}\times {\mathrm{10}}^{\mathrm{3}}]$, and for the stiffness matrix they are $\text{diag}\left({\mathit{K}}_{\mathrm{e}}\right)=[\mathrm{1.7}\times {\mathrm{10}}^{\mathrm{4}},\mathrm{6.7}\times {\mathrm{10}}^{\mathrm{4}},\mathrm{8.7}\times {\mathrm{10}}^{\mathrm{4}}]$, computed according to Eqs. (27) and (28). The coefficients k_{Ω} of each degree of freedom are obtained as ${\mathit{k}}_{\mathrm{\Omega}}=[\mathrm{1.7},\phantom{\rule{0.33em}{0ex}}\mathrm{1.4},\phantom{\rule{0.33em}{0ex}}\mathrm{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_{ω}.
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 ${\mathbf{\Phi}}_{\mathrm{1}}={\mathrm{\Phi}}_{\mathrm{1},\mathrm{x}}{\mathit{e}}_{\mathrm{x}}+{\mathrm{\Phi}}_{\mathrm{1},y}{\mathit{e}}_{\mathrm{y}}$, ${\mathbf{\Phi}}_{\mathrm{2}}={\mathrm{\Phi}}_{\mathrm{2},\mathrm{x}}{\mathit{e}}_{\mathrm{x}}+{\mathrm{\Phi}}_{\mathrm{2},y}{\mathit{e}}_{\mathrm{y}}$, and $\mathit{u}={\mathbf{\Phi}}_{\mathrm{1}}{q}_{\mathrm{1}}+{\mathbf{\Phi}}_{\mathrm{2}}{q}_{\mathrm{2}}$. The developments of appendices A and B should be used for general cases.
4.4 Twodegreesoffreedom model of a landbased or fixedbottom turbine
We consider a system of three bodies: tower (or support structure), nacelle, and rotor. The system represents a landbased wind turbine or a fixedbottom 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 function^{2} in the fore–aft direction, denoted Φ_{1}=Φ_{1}e_{x}. For hand calculations and as a first approximation, the first mode shape of a massless beam with a top mass may be used: ${\mathrm{\Phi}}_{\mathrm{1}}\left(z\right)=\mathrm{1}\mathrm{cos}(z\mathit{\pi}/L/\mathrm{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 $\mathit{q}=(q,\mathit{\psi})$, 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 towertop point remains along the x axis, neglecting the socalled nonlinear geometric effect. However, nonlinear geometric effects can be included using geometric stiffening corrections (see Appendix C or Branlard, 2019). The aerodynamic thrust and torque are denoted f_{a} and τ_{a}, respectively, and act at the rotor center (point R). The lowspeed shaft generator torque is written as τ_{g}. The distributed loads on the tower, p_{x} (from aerodynamics and hydrodynamics), are projected against the shape function to obtain the generalized forces ${f}_{\mathrm{e}}={\int}_{\mathrm{0}}^{{L}_{\mathrm{T}}}{p}_{\mathrm{x}}(z,t){\mathrm{\Phi}}_{\mathrm{1}}\left(z\right)\mathrm{d}z$. The moments of inertia of the rotor in its coordinates are $({J}_{\text{x,R}},{J}_{\oplus ,\mathrm{R}},{J}_{\oplus ,\mathrm{R}})$. We note that M_{e},K_{e}, and D_{e} are the generalized mass, stiffness, and damping, respectively, associated with a given shape function ${M}_{\mathrm{e}}={\int}_{\mathrm{0}}^{{L}_{\mathrm{T}}}m\left(z\right){\mathrm{\Phi}}_{\mathrm{1}}^{\mathrm{2}}\left(z\right)\mathrm{d}z$, ${K}_{\mathrm{e}}={\int}_{\mathrm{0}}^{{L}_{\mathrm{T}}}\text{EI}\left(z\right){\left[\frac{{d}^{\mathrm{2}}{\mathrm{\Phi}}_{\mathrm{1}}}{\mathrm{d}{z}^{\mathrm{2}}}\left(z\right)\right]}^{\mathrm{2}}\mathrm{d}z$, D_{e}=2ζM_{e}ω_{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 towertop mass (K_{gt}) and its own weight (K_{gw}) is obtained using Eq. (B5), as ${K}_{\mathrm{g}}={K}_{\text{gt}}+{K}_{\text{gw}}$, 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 (fixedbottom foundation or floater), additional geometric stiffening terms would be present (see Appendix C). The shape function frequency is obtained as
The application of the symbolic framework leads to the following equations of motion (rearranged for interpretability):
where
and
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ω^{2}r with $\mathit{\omega}={\mathit{\nu}}_{\mathrm{y}}\dot{q}$); 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 (M_{N}x_{NG}+M_{R}x_{NR}) and (M_{N}z_{NG}+M_{R}z_{NR})).
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:
For the NREL 5 MW reference turbine (Jonkman et al., 2009), the different numerical values are g=9.807 m s^{−2}, θ_{t}=5^{∘}, ${x}_{\text{NR}}=\mathrm{5.0}$ m, z_{NR}=2.4 m, L_{T}=87.6 m, z_{NG}=1.75 m, x_{NG}=1.9 m, ${M}_{\mathrm{R}}=\mathrm{1.1}\times {\mathrm{10}}^{\mathrm{5}}$ kg, ${J}_{\text{x,R}}=\mathrm{3.86}\times {\mathrm{10}}^{\mathrm{7}}$ kg m^{2}, ${J}_{\oplus ,\mathrm{R}}=\mathrm{1.92}\times {\mathrm{10}}^{\mathrm{7}}$ kg m^{2}, ${M}_{\mathrm{N}}=\mathrm{2.4}\times {\mathrm{10}}^{\mathrm{5}}$ kg, ${J}_{\text{y,N}}=\mathrm{1.01}\times {\mathrm{10}}^{\mathrm{6}}$ kg m^{2}, and ${M}_{\text{RNA}}=\mathrm{3.5}\times {\mathrm{10}}^{\mathrm{5}}$ kg. The first fore–aft shape function of the NREL 5 MW turbine tower and its derivatives are
with $\stackrel{\mathrm{\u203e}}{z}=z/L$, a_{2}=0.7004, a_{3}=2.1963, ${a}_{\mathrm{4}}=\mathrm{5.6202}$, a_{5}=6.2275, and ${a}_{\mathrm{6}}=\mathrm{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, ${M}_{\mathrm{e}}=\mathrm{5.4}\times {\mathrm{10}}^{\mathrm{4}}$, ${K}_{\mathrm{e}}=\mathrm{1.91}\times {\mathrm{10}}^{\mathrm{6}}$, ${K}_{\mathrm{g}}=\mathrm{5.2}\times {\mathrm{10}}^{\mathrm{4}}\mathrm{1.0}\times {\mathrm{10}}^{\mathrm{4}}=\mathrm{6.20}\times {\mathrm{10}}^{\mathrm{4}}$, ${\mathit{\omega}}_{\mathrm{e}}=\sqrt{({K}_{\mathrm{e}}+{K}_{\mathrm{g}})/{M}_{\mathrm{e}}}=\mathrm{5.85}$ rad s^{−1}. These numerical values, with q=0, lead to ${M}_{\mathrm{q}}=\mathrm{4.375}\times {\mathrm{10}}^{\mathrm{5}}$ and ${K}_{\mathrm{q}}=\mathrm{1.849}\times {\mathrm{10}}^{\mathrm{9}}$. The first fore–aft mode of the wind turbine has a natural frequency of $f=\sqrt{{K}_{\mathrm{q}}/{M}_{\mathrm{q}}}=\mathrm{0.3272}$ Hz. 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 Threedegreesoffreedom model of a landbased or fixedbottom 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, Φ_{1}=Φ_{1}e_{x} and Φ_{2}=Φ_{2}e_{y}. The degrees of freedom are $\mathit{q}=({q}_{\mathrm{1}},{q}_{\mathrm{2}},\mathit{\psi})$, where q_{1} and q_{2} 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.
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 f_{a} and τ_{a}, acting at point R. The distributed loads on the tower, p_{x} and p_{y} (from aerodynamics and hydrodynamics), are projected against the shape functions to obtain the generalized forces ${f}_{\mathrm{e}\mathrm{1}}=\int {\mathrm{\Phi}}_{\mathrm{1}}{p}_{\mathrm{x}}\mathrm{d}z$ and ${f}_{\mathrm{e}\mathrm{2}}=\int {\mathrm{\Phi}}_{\mathrm{2}}{p}_{\mathrm{y}}\mathrm{d}z$. The moments of inertia of the rotor in its coordinates are $({J}_{\text{x,R}},{J}_{\oplus ,\mathrm{R}},{J}_{\oplus ,\mathrm{R}})$. We note that M_{e}, K_{e}, and D_{e} are the generalized mass, stiffness, and damping, respectively, associated with a given shape function (e.g., ${M}_{\mathrm{e}\mathrm{11}}=\int {\mathrm{\Phi}}_{\mathrm{1}}^{\mathrm{2}}m\left(z\right)\mathrm{d}z$, 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 firstorder smallangle approximation for θ_{t} and a secondorder approximation for ν_{x} and ν_{y}. It is observed from Eq. (E14) that a firstorder 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 fourthorder 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 1 m in the x and y directions, and the rotational speed is 5 rpm. We report the mean relative error, ϵ, and the coefficient of determination, R^{2}, 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 smallangle approximation, which is different in OpenFAST (using the closest orthonormal matrix; Jonkman, 2009) 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.
4.6 Threedegreesoffreedom 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 $\mathit{q}=(x,\mathit{\varphi},{q}_{\mathrm{T}})$, where x is the floater surge, ϕ is the floater pitch, and q_{T} 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 p_{x}. The equations of motion are given in Appendix E3. The equations were simplified using a firstorder smallangle approximation of θ_{t} and ϕ_{y} and a secondorder approximation for ν_{y}.
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=0 m, ϕ=2^{∘}, and q_{T}=1 m. The results are presented in Fig. 8. We observe again that the results from the two models correlate to a high degree.
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.799 Hz, respectively, whereas OpenFAST returned 0.095 and 0.795 Hz. The 4 % error in the pitch frequency appears reasonable in view of the approximations used.
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 verticalaxis turbine concepts. Dedicated code can be generated for specific applications for increased performance. For instance, implicit integrators with iterative Newton–Raphsonlike 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 codesign (e.g., using methods such as linearmatrixinequalitybased 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 handbased 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 applicationspecific (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 secondorder 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 (Branlard, 2019).
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 loworder polynomial (usually linear or quadratic) approximations to eliminate highorder terms (see, e.g., Table 1 of Wallrapp, 1994). 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 “midfidelity” 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 highfidelity reducedorder 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).
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 closedloop systems and arbitrary constraints.
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 n_{e}, the flexible degrees of freedom are q_{e}, and the shape functions are gathered into a matrix Φ of size (3×n_{e}). The equations of motion, given in Eq. (12), are repeated below:
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 $\left[\stackrel{\mathrm{\u0303}}{\phantom{\rule{0.33em}{0ex}}\phantom{\rule{0.33em}{0ex}}}\right]$ represents the skew symmetric matrix. M is the mass of the body. The vector s_{CM} is the vector from the origin of the body to undeflected center or mass (CM) of the body. The notations C_{t} (n_{e}×3) and C_{r} (n_{e}×3) are introduced to match Wallrapp's notations. The vector s_{P} 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 ${\mathit{s}}_{{P}_{\mathrm{0}}}$ and the displacement field u, such that ${\mathit{s}}_{\mathrm{P}}={\mathit{s}}_{{P}_{\mathrm{0}}}+\mathit{u}$. Typically, the displacement field is given by u=Φq_{e}, but a higherorder expansion can also be introduced (see Wallrapp, 1994, and Appendix D4). Wallrapp also includes the elementary mass moment of inertia, which results in additional terms in the integrals (see Wallrapp, 1994). Such contributions are relevant, for instance, when considering the torsion of a beam (see Branlard, 2019). The block matrices M_{xx}, M_{xe}, and M_{ee} do not depend on the deformation of the body and are therefore constant. The other terms are functions of q_{e}. They may be expressed as linear combinations of constant integrals (see Appendix D3).
The quadratic velocity terms, k_{ω}, are given as
where
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 O_{e}Ω, where
The body elastic forces are given by
where K_{e} and D_{e} 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:
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, EI_{x} and EI_{y}. We assume that the displacement field is such that the shape functions are functions of z only: $\mathit{u}(z,t)={\sum}_{i=\mathrm{1}}^{{n}_{\mathrm{e}}}{\mathbf{\Phi}}_{i}\left(z\right){q}_{\mathrm{e},\mathrm{i}}\left(t\right)$. We also assume that the shape functions satisfy at least the geometric boundary conditions. The kinetic energy of the beam is $T=\frac{\mathrm{1}}{\mathrm{2}}{\int}_{\mathrm{0}}^{L}m{\dot{u}}^{\mathrm{2}}\mathrm{d}z=\frac{\mathrm{1}}{\mathrm{2}}{\sum}_{i}{\sum}_{j}{M}_{\mathrm{e},ij}{\dot{q}}_{\mathrm{e},j}{\dot{q}}_{\mathrm{e},\mathrm{i}}$, where M_{e,ij} is (see Eq. A7)
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 Branlard, 2019). The potential energy (strain energy) of the beam is obtained as $V=\frac{\mathrm{1}}{\mathrm{2}}{\sum}_{i}{\sum}_{j}{K}_{\mathrm{e},ij}{q}_{\mathrm{e},\mathrm{i}}{q}_{\mathrm{e},j}$, where K_{e,ij} represents the elements of the stiffness matrix, which, under the assumption of small deformations, are given by
Elongation and torsional strains (EA and GK_{t}) 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 δq_{e,i} provides the generalized force as (see Eq. A17)
The equations of motion of the isolated beam are then written in matrix form as
where ${\mathit{q}}_{\mathrm{e}}=[{q}_{\mathrm{e},\mathrm{1}},\mathrm{\dots},{q}_{\mathrm{e},n}]$. Damping is typically added a posteriori to the equations, where the Rayleigh damping assumption is often used: ${\mathit{D}}_{\mathrm{e}}=\mathit{\alpha}{\mathit{M}}_{\mathrm{e}}+\mathit{\beta}{\mathit{K}}_{\mathrm{e}}$ (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 ${K}_{\mathrm{e},ii}={\mathit{\omega}}_{\mathrm{e},\mathrm{i}}^{\mathrm{2}}{M}_{\mathrm{e},ii}$, with ${\mathit{\omega}}_{\mathrm{e},\mathrm{i}}=\sqrt{{K}_{\mathrm{e},ii}/{M}_{\mathrm{e},ii}}$ the eigenfrequency of the beam mode i. The modal damping is then given by ${D}_{\mathrm{e},ii}=\mathrm{2}{\mathit{\zeta}}_{\mathrm{i}}{M}_{\mathrm{e},ii}{\mathit{\omega}}_{\mathrm{e},\mathrm{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 $\mathit{n}=\frac{\partial}{\partial z}\left[N\left(z\right)\frac{\partial \mathit{u}}{\partial z}\right]$, with components in the y and z directions (see Branlard, 2019). The generalized force associated with this loading is then ${Q}_{\mathrm{N},\mathrm{i}}={\int}_{\mathrm{0}}^{L}{\mathbf{\Phi}}_{i}\cdot \mathit{n}\phantom{\rule{0.125em}{0ex}}\mathrm{d}z$. Inserting the expression of n and u, the generalized force has the form of a stiffness term: ${Q}_{\mathrm{N},\mathrm{i}}={\sum}_{j}{K}_{\mathrm{N},ij}{q}_{\mathrm{e},j}$, with
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.
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
where K_{g} is the geometric stiffness matrix of shape n_{e}×n_{e}. In general, this matrix is timedependent, 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 ($\dot{\mathit{\omega}}$), 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 (F^{k}), 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 ${\mathit{K}}_{{\mathrm{g}}^{\ast}}$, as follows:
where the indices α and β run on the x, y, and z coordinates of the body reference frame. The matrices ${\mathit{K}}_{{\mathrm{g}}^{\ast},\mathit{\alpha}}$ or ${\mathit{K}}_{{\mathrm{g}}^{\ast},\mathit{\alpha}\mathit{\beta}}$ have the shape n_{e}×n_{e} 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, K_{gt,z} is the geometric stiffness matrix corresponding to a unit acceleration in the z direction, ${\mathit{K}}_{\mathrm{g}\mathit{\omega},\text{xy}}^{k}$ is the geometric stiffness matrix corresponding to a unit gyration about the x and y directions (centrifugal effect), and ${\mathit{K}}_{\text{gF},\mathrm{x}}^{k}$ 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, $\mathrm{2}m\left(z\right)\mathit{\omega}\times \left({\sum}_{j}{\mathbf{\Phi}}_{j}\right(z\left){\dot{q}}_{\mathrm{e},j}\right)$, is proportional to ${\dot{\mathit{q}}}_{\mathrm{e}}$. Because the instantaneous beam slope is proportional to q_{e}, the geometric stiffening term consists of terms of the form ${q}_{\mathrm{e},j}{\dot{q}}_{\mathrm{e},k}$. In Table 1 of Wallrapp (1994), it is stated that nonlinear terms of the form ${q}_{\mathrm{e},j}{\dot{q}}_{\mathrm{e},k}$ are neglected. Yet, if the steadystate deflection q_{e} 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 q_{e} and ${\dot{\mathit{q}}}_{\mathrm{e}}$ (expressions are provided in Eq. C2). We note that the terms ${\mathit{K}}_{{\mathrm{g}}^{\ast}}$ have different units; for instance, the terms ${\mathit{K}}_{\text{gt}{,}^{\ast}}$ are expressed in N s^{2} 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 finiteelement 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:
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 K_{gt,z} is obtained by considering a unit vertical acceleration:
We write z_{k} as the coordinate of node k along the beam. The ij component of the matrix ${\mathit{K}}_{\text{gF},z}^{k}$ is obtained as
The ij component of the matrix K_{gω,αβ} is obtained by considering unit centrifugal loads generated using independent rotations around the unit vectors e_{x}, e_{y}, and e_{z}:
Similarly, the ij component of the matrix K_{gr,α} is
The Coriolis force, $\mathrm{2}m\left(z\right)\mathit{\omega}\times \left({\sum}_{j}{\mathbf{\Phi}}_{j}\right(z\left){\dot{q}}_{\mathrm{e},j}\right)$, can also have an axial contribution. Using Eq. (C3), the generalized force is
Equation (C8) may be rearranged by introducing a threedimensional 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 ${\mathit{k}}_{\mathit{\sigma},\text{Cor}}({\mathit{q}}_{\mathrm{e}},{\dot{\mathit{q}}}_{\mathrm{e}})$ is added to the righthand side of Eq. (C1).
C3 Integration into the equations of motion
The term k_{σ}=K_{g}q_{e} 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 (${\mathit{M}}_{{\mathrm{e}}^{\ast}}$), the term k_{ω,e}, and the term f_{e}, respectively. For instance, the term $\sum {a}_{\mathit{\alpha}}{\mathit{K}}_{\text{gt},\mathit{\alpha}}{\mathit{q}}_{\mathrm{e}}$ can be reorganized as [K_{gt}]q_{e}⋅a (using loose notations); therefore, the mass matrix can be updated such that M_{xe} becomes M_{xe}+[K_{gt}]q_{e}. When a Taylor expansion is used, such integration is easily implemented as a firstorder term (see Appendix D3).
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., Lemmer, 2018, and Branlard, 2019), the Jacobians are sometimes stacked into a matrix form:
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 treelike 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 timevarying 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 (Wallrapp, 1994; Schwertassek and Wallrapp, 1999; Lemmer, 2018). In such a formulation, the linear rotation field is obtained as $\mathit{I}+\stackrel{\mathrm{\u0303}}{\mathbf{\Psi}\mathit{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, ${\mathit{s}}_{\mathrm{P}}={\mathit{s}}_{{P}_{\mathrm{0}}}+\mathit{u}$, where the displacement field is u=Φq_{e}. The undeflected position of the structure (${\mathit{s}}_{{P}_{\mathrm{0}}}$) is constant, and the shape functions are known at the initialization; the only timevarying terms are the degrees of freedom q_{e}. 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 q_{e}. The precomputed integrals are referred to as “shape integrals”. For a given term T (standing, for instance, for M_{θ,θ}, C_{t}, C_{r}, G_{r}, G_{e}, or O_{e}), the shape integral expansion is
If T is an array, T^{0} and ${\mathit{T}}_{j}^{\mathrm{1}}$ have the same shape as T. As an example, the application of the shape integral expansion to the term M_{xθ} (see Eq. A3) gives
with
The zeroth and firstorder shape integrals always consist of integrals over the components of ${\mathit{s}}_{{P}_{\mathrm{0}}}$ 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 Wallrapp, 1994, or the σ, Σ, Υ, and Ψ terms introduced by Shabana, 2013). The zeroth and firstorder shape integrals are stored using a Taylor objectoriented 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 finiteelement formulation (see Schwertassek and Wallrapp, 1999).
The geometric stiffness introduced in Appendix C is linear in the elastic degrees of freedom q_{e}. Therefore, the unit geometric stiffness matrices (which are also shape integrals) can be conveniently added into the firstorder terms of Eq. (D2). For instance, if we write M_{ex} (given in Eq. A6) using a firstorder expansion, ${\mathit{M}}_{\text{ex}}={\mathit{M}}_{\text{ex}}^{\mathrm{0}}+{\mathit{M}}_{\text{ex}}^{\mathrm{1}}{\mathit{q}}_{\mathrm{e}}$, then the geometric stiffening effect can directly be inserted into the firstorder term, such that ${\mathit{M}}_{\text{ex}}^{\mathrm{1}}$ becomes ${\mathit{M}}_{\text{ex}}^{\mathrm{1}}+{\mathit{K}}_{\text{gt}}$. Similarly, the term K_{gr} can be inserted into ${\mathit{M}}_{\mathit{\theta}\mathrm{e}}^{\mathrm{1}}$, K_{gω} into ${\mathit{O}}_{\mathrm{e}}^{\mathrm{1}}$, K_{gF} into Φ^{1}, and K_{gτ} 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 firstorder coefficient is that it could make the mass matrix symmetric (if the user code assumes ${\mathit{M}}_{\text{xe}}={\mathit{M}}_{\text{ex}}^{t}$), 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 (Wallrapp, 1993, 1994), the displacement field is assumed to be a function of the degrees of freedom, u=Φ_{u}(q_{e})q_{e}, 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 socalled “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 treelike 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:
where ${\mathit{J}}_{{v}_{\mathrm{P}},r}$ and ${\mathit{J}}_{{\mathit{\omega}}_{\mathrm{P}},r}$ are the partial velocities of point P with respect to the degree of freedom r, f_{P} 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 h_{r} 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 f_{P} 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”):
The terms f_{P,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 righthand side of the equation of motion (Eq. 23) are obtained as
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.
In this section, we present the equations of motion for the examples presented in Sect. 4.
E1 Twodegreesoffreedom model of a landbased or fixedbottom 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, $v{\widehat{t}}_{\mathrm{x}}$ is the unit vector in the x direction of the tower frame. The degrees of freedom are $\mathit{q}=(q,\mathit{\psi})$. The kinematics of the tower (at its origin) are zero:
All Jacobians are zero except ${\mathit{J}}_{\mathrm{e},\mathrm{1}T}=\mathrm{1}$. The inertial force, torque, and elastic force are
The nacelle kinematics (at its center of mass) are
The Jacobians with respect to q are
The inertial force and torque on the nacelle are
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 Threedegreesoffreedom model of a landbased or fixedbottom wind turbine
The equations of motion for the model presented in Sect. 4.5, with $\mathit{q}=({q}_{\mathrm{1}},{q}_{\mathrm{2}},\mathit{\psi})$, are given in this section. The elements of the mass matrix are
The elements of the forcing vector are
E3 Threedegreesoffreedom model of a floating wind turbine
The equations of motion for the model presented in Sect. 4.6, with $\mathit{q}=(x,\mathit{\varphi},{q}_{\mathrm{T}})$, are given in this section. The elements of the mass matrix are
The elements of the forcing vector are
The latest source code of YAMS is available on GitHub as a subpackage of the Wind Energy Library, WELIB (http://github.com/ebranlard/welib/, Branlard, 2022a). A static version is available on Zenodo (https://doi.org/10.5281/zenodo.7306075, Branlard, 2022b). The examples given in this articles are found in the folder welib/yams/papers of the repository.
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.
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.
This work was funded under the Technology Commercialization Fund Project, supported by the DOE's Wind Energy Technologies Office.
This paper was edited by Raimund Rolfes and reviewed by two anonymous referees.
ANSYS: https://www.ansys.com/ (last access: 19 March 2022), 2022. a
Bauchau, O. A.: Flexible Multibody Dynamics, Solid Mechanics and Its Applications, Springer, Dordrecht, https://doi.org/10.1007/9789400703353, 2011. a
Benner, P., Gugercin, S., and Willcox, K.: A Survey of ProjectionBased Model Reduction Methods for Parametric Dynamical Systems, SIAM Rev., 57, 483–531, https://doi.org/10.1137/130932715, 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 RayleighRitz approximation: The general framework behind and beyond Flex, Wind Energy, 22, 877–893, https://doi.org/10.1002/we.2327, 2019. a, b, c, d, e, f, g, h, i, j, k, l, m
Branlard, E.: WELIB, Wind Energy Library, GitHub repository http://github.com/ebranlard/welib/ (last access: November 2022), 2022a. a, b, c
Branlard, E.: WELIB, Zenodo, https://doi.org/10.5281/zenodo.7306075, 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 landbased wind turbine: a step towards digitaltwin simulations, Wind Energ. Sci., 5, 1155–1167, https://doi.org/10.5194/wes511552020, 2020a. a
Branlard, E., Jonkman, J., Dana, S., and Doubrawa, P.: A digital twin based on OpenFAST linearizations for realtime load and fatigue estimation of landbased turbines, J. Phys. Conf. Ser., 1618, 022030, https://doi.org/10.1088/17426596/1618/2/022030, 2020b. a
Docquier, N., Poncelet, A., and Fisette, P.: ROBOTRAN: a powerful symbolic gnerator of multibody models, Mech. Sci., 4, 199–219, https://doi.org/10.5194/ms41992013, 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, https://doi.org/10.1115/DETC201313470, 2013. a, b, c
Geisler, J.: CADynTub: Wind Turbine Model from OpenFAST Data using CADyn Equations of Motion, Github, https://github.com/jgeisler0303/CADynTurb (last access: November 2022), 2021. a
Géradin, M. and Cardona, A.: Flexible Multibody Dynamics: A Finite Element Approach, Wiley, New York, ISBN 9780471489900, 2001. a
Jelenić, G. and Crisfield, M.: Geometrically exact 3D beam theory: implementation of a straininvariant finite element for statics and dynamics, Comput. Method. Appl. M., 171, 141–171, https://doi.org/10.1016/S00457825(98)002497, 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. Opensource wind turbine simulation tool, Zenodo, https://doi.org/10.5281/zenodo.6324288, 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/TP50038060, National Renewable Energy Laboratory, https://doi.org/10.2172/947422, 2009. a, b
Jonkman, J. M.: Dynamics of offshore floating wind turbines–model development and verification, Wind Energy, 12, 459–492, https://doi.org/10.1002/we.347, 2009. a
Jonkman, J. M., Branlard, E. S. P., and Jasa, J. P.: Influence of wind turbine design parameters on linearized physicsbased models in OpenFAST, Wind Energ. Sci., 7, 559–571, https://doi.org/10.5194/wes75592022, 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, https://doi.org/10.1137/0113030, 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 1619, 2009. a, b
Kurtz, T., Eberhard, P., Henninger, C., and Schiehlen, W.: From Neweul to NeweulM2: symbolical equations of motion for multibody system analysis and synthesis, Multibody Syst.Dyn., 24, 25–41, https://doi.org/10.1007/s110440109187x, 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, SaintHubert, Quebéc, 31 May–1 June 2007. a
Lemmer, F.: Loworder modeling, controller design and optimization of floating offshore wind turbines, PhD thesis, Universität Stuttgart, http://elib.unistuttgart.de/handle/11682/10543 (last access: November 2022), 2018. a, b, c, d
MBDyn: https://www.mbdyn.org/ (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, http://www.motiongenesis.com (last access: November 2022), 2016. a
Øye, S.: Fix Dynamisk, aeroelastisk beregning af vindmøllevinger, Report AFM8308, 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 inequalitybased control design, Wind Energy, 23, 1792–1809, https://doi.org/10.1002/we.2516, 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, https://doi.org/10.1088/17426596/659/1/012001, 2015. a
Simo, J.: A finite strain beam formulation. The threedimensional dynamic problem. Part I, Comput. Method. Appl. M., 49, 55–70, https://doi.org/10.1016/00457825(85)900507, 1985. a
SIMPACK: https://www.3ds.com/productsservices/simulia/products/simpack/ (last access: 19 March 2022), 2022. a
Sønderby, I. and Hansen, M. H.: Openloop frequency response analysis of a wind turbine using a highorder linear aeroelastic model, Wind Energy, 17, 1147–1167, https://doi.org/10.1002/we.1624, 2014. a
Steindl, A. and Troger, H.: Methods for dimension reduction and their application in nonlinear dynamics, Int. J. Solids Struct., 38, 2131–2147, https://doi.org/10.1016/S00207683(00)001578, 2001. a
SymPy: https://www.sympy.org (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, https://doi.org/10.1007/s11071021066939, 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, https://doi.org/10.1007/9789401706254_33, 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, https://doi.org/10.1002/nme.1620320818, 10.1002/(ISSN)10970207, 1991. a
 Abstract
 Copyright statement
 Introduction
 Method to obtain the equations of motion
 Implementation into a symbolic framework
 Wind energy applications
 Discussions
 Conclusions
 Appendix A: Equations for a flexible body and shape integrals
 Appendix B: Application of the shape function approach to an isolated beam
 Appendix C: Geometric stiffness
 Appendix D: Alternative formulations
 Appendix E: Equations of motion of simple wind turbine models
 Code availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References
 Abstract
 Copyright statement
 Introduction
 Method to obtain the equations of motion
 Implementation into a symbolic framework
 Wind energy applications
 Discussions
 Conclusions
 Appendix A: Equations for a flexible body and shape integrals
 Appendix B: Application of the shape function approach to an isolated beam
 Appendix C: Geometric stiffness
 Appendix D: Alternative formulations
 Appendix E: Equations of motion of simple wind turbine models
 Code availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References