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

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 (Ø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 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 Hansen, 2014) 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 (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 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 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 co-design approaches. We implemented an open-source 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 closed-source 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 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 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(\dot{}\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)

$$\begin{array}{}\text{(1)}& {\displaystyle }{\mathit{r}}_{\mathrm{P}}={\mathit{r}}_{\mathrm{i}}+{\mathit{s}}_{\mathrm{P}}={\mathit{r}}_{\mathrm{i}}+{\mathit{s}}_{P_{\mathrm{0}}}+{\mathit{u}}_{\mathrm{P}},\text{(2)}& {\displaystyle }{\mathit{v}}_{\mathrm{P}}={\mathit{v}}_{\mathrm{i}}+{\mathit{\omega }}_{\mathrm{i}}\times {\mathit{s}}_{\mathrm{P}}+({\dot{\mathit{u}}}_{\mathrm{P}}{)}_{\mathrm{i}},\text{(3)}& {\displaystyle }\begin{array}{rl}{\mathit{a}}_{\mathrm{P}}=& {\mathit{a}}_{\mathrm{i}}+{\mathit{\omega }}_{\mathrm{i}}\times ({\mathit{\omega }}_{\mathrm{i}}\times {\mathit{s}}_{\mathrm{P}})+{\dot{\mathit{\omega }}}_{\mathrm{i}}\times {\mathit{s}}_{\mathrm{P}}\\ & +\mathrm{2}{\mathit{\omega }}_{\mathrm{i}}\times ({\dot{\mathit{u}}}_{\mathrm{P}}{)}_{\mathrm{i}}+({\ddot{\mathit{u}}}_{\mathrm{P}}{)}_{\mathrm{i}},\end{array}\end{array}$$

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₀ 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{}\right)$ and $(\dot{}{)}_{\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{̃}}{\mathit{\omega }}}_{\mathrm{i}}\right]={\dot{\mathit{R}}}_{\mathrm{i}}{\mathit{R}}_{\mathrm{i}}^{\mathrm{T}}$, where $\left[\stackrel{\mathrm{̃}}{}\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 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 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:

$$\begin{array}{}\text{(4)}& f_r+f_r^{\ast }=\mathrm{0},r=\mathrm{1}\mathrm{\dots }{n}_{\mathrm{q}},\end{array}$$

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:

$$\begin{array}{}\text{(5)}& f_r=\sum _{i=\mathrm{1}}^{n_{\mathrm{b}}}{f}_{\text{ri}},f_r^{\ast }=\sum _{i=\mathrm{1}}^{n_{\mathrm{b}}}{f}_{\text{ri}}^{\ast }.\end{array}$$

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

$$\begin{array}{}\text{(6)}& {\mathit{f}}_{\mathrm{i}}^{\ast }=-m_{\mathrm{i}}{\mathit{a}}_{\mathrm{G},\mathrm{i}},{\mathit{\tau }}_{\mathrm{i}}^{\ast }=-{\mathit{I}}_{\mathrm{G},\mathrm{i}}\cdot {\dot{\mathit{\omega }}}_{\mathrm{i}}-{\mathit{\omega }}_{\mathrm{i}}\times ({\mathit{I}}_{\mathrm{G},\mathrm{i}}\cdot {\mathit{\omega }}_{\mathrm{i}}),\end{array}$$

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

$$\begin{array}{}\text{(7)}& f_{\text{ri}}^{\ast }={\mathit{J}}_{\mathrm{v},\text{ri}}\cdot {\mathit{f}}_{\mathrm{i}}^{\ast }+{\mathit{J}}_{\mathit{\omega },\text{ri}}\cdot {\mathit{\tau }}_{\mathrm{i}}^{\ast },\end{array}$$

with

$$\begin{array}{}\text{(8)}& {\mathit{J}}_{\mathrm{v},\text{ri}}={\displaystyle \frac{\partial {\mathit{v}}_{\mathrm{G},\mathrm{i}}}{\partial {\dot{q}}_{\mathrm{r}}}},{\mathit{J}}_{\mathit{\omega },\text{ri}}={\displaystyle \frac{\partial {\mathit{\omega }}_{\mathrm{i}}}{\partial {\dot{q}}_{\mathrm{r}}}},\end{array}$$

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

$$\begin{array}{}\text{(9)}& f_{\text{ri}}={\mathit{J}}_{\mathrm{v},\text{ri}}\cdot {\mathit{f}}_{\mathrm{i}}+{\mathit{J}}_{\mathit{\omega },\text{ri}}\cdot {\mathit{\tau }}_{\mathrm{i}}.\end{array}$$

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:

$$\begin{array}{}\text{(10)}& f_{\text{ri}}=\sum _j{\displaystyle \frac{\partial {\mathit{v}}_{P_j}}{\partial {\dot{q}}_{\mathrm{r}}}}\cdot {\mathit{f}}_{i,j}+\sum _k{\mathit{J}}_{\mathit{\omega },\text{ri}}\cdot {\mathit{\tau }}_{i,k},\end{array}$$

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 zeroth-order¹ Rayleigh–Ritz approximation, the displacement field at a given point, P, is given by the sum of shape function contributions:

$$\begin{array}{}\text{(11)}& \mathit{u}\left(P\right)=\sum _{j=\mathrm{1}}^{n_{\mathrm{e},i}}{\mathbf{\Phi }}_{ij}\left(P\right){\mathit{q}}_{\mathrm{e},ij}\left(t\right),\end{array}$$

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)

$$\begin{array}{}\text{(12)}& {\displaystyle }\begin{array}{rl}{\left[\begin{array}{ccc}{\mathit{M}}_{\text{xx}}^{\prime }& {\mathit{M}}_{\text{x}\mathit{\theta }}^{\prime }& {\mathit{M}}_{\text{xe}}^{\prime }\\ & {\mathit{M}}_{\mathit{\theta }\mathit{\theta }}^{\prime }& {\mathit{M}}_{\mathit{\theta }\text{e}}^{\prime }\\ \text{sym.}& & {\mathit{M}}_{\text{ee}}^{\prime }\end{array}\right]}_{\mathrm{i}}\left[\begin{array}{c}{\mathit{a}}_{\mathrm{i}}^{\prime }\\ {\dot{\mathit{\omega }}}_{\mathrm{i}}^{\prime }\\ {\ddot{\mathit{q}}}_{\mathrm{e},\mathrm{i}}\end{array}\right]& +{\left[\begin{array}{c}{\mathit{k}}_{\mathit{\omega },\text{x}}^{\prime }\\ {\mathit{k}}_{\mathit{\omega },\mathit{\theta }}^{\prime }\\ {\mathit{k}}_{\mathit{\omega },\text{e}}^{\prime }\end{array}\right]}_{\mathrm{i}}\\ & +{\left[\begin{array}{c}\mathrm{0}\\ \mathrm{0}\\ {\mathit{k}}_{\mathrm{e}}\end{array}\right]}_{\mathrm{i}}={\left[\begin{array}{c}{\mathit{f}}_{\mathrm{x}}^{\prime }\\ {\mathit{f}}_{\mathit{\theta }}^{\prime }\\ {\mathit{f}}_{\mathrm{e}}\end{array}\right]}_{\mathrm{i}},\end{array}\end{array}$$

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

$$\begin{array}{}\text{(13)}& {\displaystyle }{\mathit{f}}_{\mathrm{i}}^{\ast }& {\displaystyle }=-{\mathit{R}}_{\mathrm{i}}\left[{\mathit{M}}_{\text{xx}}^{\prime }{\mathit{a}}_{\mathrm{i}}^{\prime }+{\mathit{M}}_{\mathrm{x}\mathit{\theta }}^{\prime }{\dot{\mathit{\omega }}}_{\mathrm{i}}^{\prime }+{\mathit{M}}_{\text{xe}}{\ddot{\mathit{q}}}_{\mathrm{e},\mathrm{i}}+{\mathit{k}}_{\mathit{\omega },\mathrm{x}}^{\prime }\right],\text{(14)}& {\displaystyle }{\mathit{\tau }}_{\mathrm{i}}^{\ast }& {\displaystyle }=-{\mathit{R}}_{\mathrm{i}}\left[{\mathit{M}}_{\mathit{\theta }x}^{\prime }{\mathit{a}}_{\mathrm{i}}^{\prime }+{\mathit{M}}_{\mathit{\theta }\mathit{\theta }}^{\prime }{\dot{\mathit{\omega }}}_{\mathrm{i}}^{\prime }+{\mathit{M}}_{\mathit{\theta }\mathrm{e}}{\ddot{\mathit{q}}}_{\mathrm{e},\mathrm{i}}+{\mathit{k}}_{\mathit{\omega },\mathit{\theta }}^{\prime }\right],\text{(15)}& {\displaystyle }{\mathit{h}}_{\mathrm{i}}^{\ast }& {\displaystyle }=-\left[{\mathit{M}}_{\text{ex}}^{\prime }{\mathit{a}}_{\mathrm{i}}^{\prime }+{\mathit{M}}_{\mathrm{e}\mathit{\theta }}^{\prime }{\dot{\mathit{\omega }}}_{\mathrm{i}}^{\prime }+{\mathit{M}}_{\text{ee}}{\ddot{\mathit{q}}}_{\mathrm{e},\mathrm{i}}+{\mathit{k}}_{\mathit{\omega },\mathrm{e}}^{\prime }\right].\end{array}$$

The external and elastic loads are

$$\begin{array}{}\text{(16)}& {\displaystyle }{\mathit{f}}_{\mathrm{i}}& {\displaystyle }={\mathit{R}}_{\mathrm{i}}{\mathit{f}}_{\mathrm{x}}^{\prime },\text{(17)}& {\displaystyle }{\mathit{\tau }}_{\mathrm{i}}& {\displaystyle }={\mathit{R}}_{\mathrm{i}}{\mathit{f}}_{\mathit{\theta }}^{\prime },\text{(18)}& {\displaystyle }{\mathit{h}}_{\mathrm{i}}& {\displaystyle }={\mathit{f}}_{\mathrm{e}}-{\mathit{k}}_{\mathrm{e}}.\end{array}$$

The components of $f_{\text{ri}}^{\ast }$ and f_ri, for r=1…n_q, are then defined as

$$\begin{array}{}\text{(19)}& {\displaystyle }& {\displaystyle }{f}_{\text{ri}}^{\ast }={\mathit{J}}_{\mathrm{v},\text{ri}}\cdot {\mathit{f}}_{\mathrm{i}}^{\ast }+{\mathit{J}}_{\mathit{\omega },\text{ri}}\cdot {\mathit{\tau }}_{\mathrm{i}}^{\ast }+{\mathit{J}}_{\mathrm{e},\text{ri}}\cdot {\mathit{h}}_{\mathrm{i}}^{\ast },\text{(20)}& {\displaystyle }& {\displaystyle }{f}_{\text{ri}}={\mathit{J}}_{\mathrm{v},\text{ri}}\cdot {\mathit{f}}_{\mathrm{i}}+{\mathit{J}}_{\mathit{\omega },\text{ri}}\cdot {\mathit{\tau }}_{\mathrm{i}}+{\mathit{J}}_{\mathrm{e},\text{ri}}\cdot {\mathit{h}}_{\mathrm{i}},\end{array}$$

with

$$\begin{array}{}\text{(21)}& {\mathit{J}}_{\mathrm{v},\text{ri}}={\displaystyle \frac{\partial {\mathit{v}}_{\mathrm{O},\mathrm{i}}}{\partial {\dot{q}}_{\mathrm{r}}}},{\mathit{J}}_{\mathit{\omega },\text{ri}}={\displaystyle \frac{\partial {\mathit{\omega }}_{\mathrm{i}}}{\partial {\dot{q}}_{\mathrm{r}}}},{\mathit{J}}_{\mathrm{e},\text{ri}}={\displaystyle \frac{\partial {\mathit{q}}_{\mathrm{e},\mathrm{i}}}{\partial q_{\mathrm{r}}}},\end{array}$$

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

$$\begin{array}{}\text{(22)}& \mathbf{e}(\mathit{q},\dot{\mathit{q}},\ddot{\mathit{q}},\mathit{u},t)=\mathit{f}+{\mathit{f}}^{\ast }=\mathit{F}(\mathit{q},\dot{\mathit{q}},\mathit{u},t)-\mathit{M}\left(\mathit{q}\right)\ddot{\mathit{q}}=\mathbf{0}\end{array}$$

or

$$\begin{array}{}\text{(23)}& \mathit{M}\left(\mathit{q}\right)\ddot{\mathit{q}}=\mathit{F}(\mathit{q},\dot{\mathit{q}},\mathit{u},t),\end{array}$$

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 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 $\dot{\mathit{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

$$\begin{array}{}\text{(24)}& \mathbf{e}({\mathit{q}}_{\mathrm{0}},{\dot{\mathit{q}}}_{\mathrm{0}},{\ddot{\mathit{q}}}_{\mathrm{0}},{\mathit{u}}_{\mathrm{0}},t)=\mathbf{0}.\end{array}$$

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

$$\begin{array}{}\text{(25)}& \begin{array}{rl}{\mathit{M}}_{\mathrm{0}}\left({\mathit{q}}_{\mathrm{0}}\right)\mathit{\delta }\ddot{\mathit{q}}& +{\mathit{C}}_{\mathrm{0}}({\mathit{q}}_{\mathrm{0}},{\dot{\mathit{q}}}_{\mathrm{0}},{\mathit{u}}_{\mathrm{0}})\mathit{\delta }\dot{\mathit{q}}\\ & +{\mathit{K}}_{\mathrm{0}}({\mathit{q}}_{\mathrm{0}},{\dot{\mathit{q}}}_{\mathrm{0}},{\ddot{\mathit{q}}}_{\mathrm{0}},{\mathit{u}}_{\mathrm{0}})\mathit{\delta }\mathit{q}={\mathit{Q}}_{\mathrm{0}}({\mathit{q}}_{\mathrm{0}},{\dot{\mathit{q}}}_{\mathrm{0}},{\mathit{u}}_{\mathrm{0}})\mathit{\delta }\mathit{u},\end{array}\end{array}$$

with

$$\begin{array}{}\text{(26)}& \begin{array}{rl}& {\mathit{M}}_{\mathrm{0}}=-{\left.{\displaystyle \frac{\partial \mathbf{e}}{\partial \ddot{\mathit{q}}}}\right|}_{\mathrm{0}},{\mathit{C}}_{\mathrm{0}}=-{\left.{\displaystyle \frac{\partial \mathbf{e}}{\partial \dot{\mathit{q}}}}\right|}_{\mathrm{0}},{\mathit{K}}_{\mathrm{0}}=-{\left.{\displaystyle \frac{\partial \mathbf{e}}{\partial \mathit{q}}}\right|}_{\mathrm{0}},\\ & {\mathit{Q}}_{\mathrm{0}}={\left.{\displaystyle \frac{\partial \mathbf{e}}{\partial \mathit{u}}}\right|}_{\mathrm{0}},\end{array}\end{array}$$

where M₀, C₀, and K₀ are the linear mass, damping, and stiffness matrices, respectively; Q₀δu is the linear forcing vector (Q₀ is the input matrix); δ indicates a small perturbation of the quantities; and |₀ 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 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 (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 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 1 Typical 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 $\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 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 (Branlard, 2022a). 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 (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 Φ₁, 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):

$$\begin{array}{}\text{(27)}& {\displaystyle }{M}_{\mathrm{e}}& {\displaystyle }=\underset{\mathrm{0}}{\overset{R}{\int }}m\left(r\right){\mathrm{\Phi }}^{\mathrm{2}}\left(r\right)\mathrm{d}r,\text{(28)}& {\displaystyle }{K}_{\mathrm{e}}& {\displaystyle }=\underset{\mathrm{0}}{\overset{R}{\int }}{\text{EI}}_{\mathrm{y}}\left(r\right){\left[{\displaystyle \frac{{\mathrm{d}}^{\mathrm{2}}\mathrm{\Phi }}{\mathrm{d}{r}^{\mathrm{2}}}}\left(r\right)\right]}^{\mathrm{2}}\mathrm{d}r.\end{array}$$

Figure 2 Sketch 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):

$$\begin{array}{}\text{(29)}& f_{\mathrm{e}}=\underset{\mathrm{0}}{\overset{R}{\int }}{p}_{\mathrm{x}}(r,t)\mathrm{\Phi }\left(r\right)\mathrm{d}r.\end{array}$$

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:

$$\begin{array}{}\text{(30)}& N\left(r\right)=\underset{\mathrm{r}}{\overset{R}{\int }}m\left(r^{\prime }\right){\mathrm{\Omega }}^{\mathrm{2}}{r}^{\prime }\mathrm{d}{r}^{\prime }.\end{array}$$

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

$$\begin{array}{}\text{(31)}& \begin{array}{rl}{K}_{\mathrm{g}}\left(\mathrm{\Omega }\right)& =\underset{\mathrm{0}}{\overset{R}{\int }}N\left(r\right){\left[{\displaystyle \frac{\mathrm{d}\mathrm{\Phi }}{\mathrm{d}r}}\right]}^{\mathrm{2}}\mathrm{d}r\\ & ={\mathrm{\Omega }}^{\mathrm{2}}\underset{\mathrm{0}}{\overset{R}{\int }}\underset{\mathrm{r}}{\overset{R}{\int }}m\left(r^{\prime }\right)r^{\prime }\mathrm{d}{r}^{\prime }{\left[{\displaystyle \frac{\mathrm{d}\mathrm{\Phi }}{\mathrm{d}r}}\right]}^{\mathrm{2}}\mathrm{d}r.\end{array}\end{array}$$

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:

$$\begin{array}{}\text{(32)}& \begin{array}{rl}{\mathit{\omega }}_{\mathrm{0}}\left(\mathrm{\Omega }\right)& =\sqrt{{\displaystyle \frac{(K_{\mathrm{e}}+K_{\mathrm{g}}(\mathrm{\Omega }\left)\right)}{M_{\mathrm{e}}}}}=\sqrt{{\mathit{\omega }}_{\mathrm{0}}^{\mathrm{2}}\left(\mathrm{0}\right)+{\displaystyle \frac{K_{\mathrm{g}}\left(\mathrm{\Omega }\right)}{M_{\mathrm{e}}}}}\\ & =\sqrt{{\mathit{\omega }}_{\mathrm{0}}^{\mathrm{2}}\left(\mathrm{0}\right)+k_{\mathrm{\Omega }}{\mathrm{\Omega }}^{\mathrm{2}}},\end{array}\end{array}$$

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 Φ₂=Φ₂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_θ=Φ₂q₂. 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

$$\begin{array}{}\text{(33)}& \begin{array}{rl}{\mathit{\omega }}_{\mathrm{0}}\left(\mathrm{\Omega }\right)& =\sqrt{{\displaystyle \frac{(K_{\mathrm{e}}+K_{\mathrm{g}}(\mathrm{\Omega })+K_{\mathit{\omega }}(\mathrm{\Omega }\left)\right)}{M_{\mathrm{e}}}}}\\ & =\sqrt{{\mathit{\omega }}_{\mathrm{0}}^{\mathrm{2}}\left(\mathrm{0}\right)+(k_{\mathrm{\Omega }}-\mathrm{1}){\mathrm{\Omega }}^{\mathrm{2}}},\end{array}\end{array}$$

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},\mathrm{1.4},\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_ω.

Figure 3 Variation 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 ${\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.

Figure 4 Model 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 function² in the fore–aft direction, denoted Φ₁=Φ₁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}-\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 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 Branlard, 2019). The aerodynamic thrust and torque are denoted f_a 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, 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 tower-top 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

$$\begin{array}{}\text{(34)}& {\displaystyle }{K}_{\text{gt}}& {\displaystyle }=-g\underset{\mathrm{0}}{\overset{L_{\mathrm{T}}}{\int }}(M_{\mathrm{R}}+M_{\mathrm{N}}){\left[{\displaystyle \frac{\mathrm{d}{\mathrm{\Phi }}_{\mathrm{1}}}{\mathrm{d}z}}\left(z\right)\right]}^{\mathrm{2}}\mathrm{d}z,\text{(35)}& {\displaystyle }{K}_{\text{gw}}& {\displaystyle }=-g\underset{\mathrm{0}}{\overset{L_{\mathrm{T}}}{\int }}\left[\underset{\mathrm{z}}{\overset{L_{\mathrm{T}}}{\int }}m\left(z^{\prime }\right)\mathrm{d}{z}^{\prime }\right]{\left[{\displaystyle \frac{\mathrm{d}{\mathrm{\Phi }}_{\mathrm{1}}}{\mathrm{d}z}}\left(z\right)\right]}^{\mathrm{2}}\mathrm{d}z.\end{array}$$

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

$$\begin{array}{}\text{(36)}& {\mathit{\omega }}_{\mathrm{e}}=\sqrt{(K_{\mathrm{e}}+K_{\mathrm{g}})/M_{\mathrm{e}}}.\end{array}$$

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

$$\begin{array}{}\text{(37)}& \left[\begin{array}{cc}{M}_{\mathrm{q}}& \mathrm{0}\\ \mathrm{0}& J_{\text{x,R}}\end{array}\right]\left[\begin{array}{c}\ddot{q}\\ \ddot{\mathit{\psi }}\end{array}\right]=\left[\begin{array}{c}{f}_{\mathrm{q}}\\ {\mathit{\tau }}_{\mathrm{a}}-{\mathit{\tau }}_{\mathrm{g}}\end{array}\right],\end{array}$$

where

$$\begin{array}{ll}\text{(38)}& {\displaystyle }{M}_{\mathrm{q}}=& {\displaystyle }{M}_{\mathrm{e}}+M_{\mathrm{N}}+M_{\mathrm{R}}{\displaystyle }& {\displaystyle }+\left(J_{\text{yN}}+J_{\oplus ,\mathrm{R}}+M_{\mathrm{N}}\left(x_{\text{NG}}^{\mathrm{2}}+z_{\text{NG}}^{\mathrm{2}}\right)\right.\\ \text{(39)}& {\displaystyle }& {\displaystyle }\left.+M_{\mathrm{R}}\left(x_{\text{NR}}^{\mathrm{2}}+z_{\text{NR}}^{\mathrm{2}}\right)\right){\mathit{\nu }}_{\mathrm{y}}^{\mathrm{2}}{\displaystyle }& {\displaystyle }+\mathrm{2}\left[(M_{\mathrm{N}}{z}_{\text{NG}}+M_{\mathrm{R}}{z}_{\text{NR}})\cos \left({\mathit{\nu }}_{\mathrm{y}}q\right)\right.\\ \text{(40)}& {\displaystyle }& {\displaystyle }\left.-(M_{\mathrm{N}}{x}_{\text{NG}}+M_{\mathrm{R}}{x}_{\text{NR}})\sin \left({\mathit{\nu }}_{\mathrm{y}}q\right)\right]{\mathit{\nu }}_{\mathrm{y}}\end{array}$$

and

$$\begin{array}{ll}\text{(41)}& {\displaystyle }{f}_{\mathrm{q}}=& {\displaystyle }{f}_{\mathrm{e}}-(K_{\mathrm{e}}+K_{\mathrm{g}})q-D_{\mathrm{e}}\dot{q}{\displaystyle }& {\displaystyle }+g{\mathit{\nu }}_{\mathrm{y}}\left[(M_{\mathrm{N}}{x}_{\text{NG}}+M_{\mathrm{R}}{x}_{\text{NR}})\cos \left({\mathit{\nu }}_{\mathrm{y}}q\right)\right.\\ \text{(42)}& {\displaystyle }& {\displaystyle }\left.+(M_{\mathrm{N}}{z}_{\text{NG}}+M_{\mathrm{R}}{z}_{\text{NR}})\sin \left({\mathit{\nu }}_{\mathrm{y}}q\right)\right]{\displaystyle }& {\displaystyle }+{\mathit{\nu }}_{\mathrm{y}}^{\mathrm{2}}{\dot{q}}^{\mathrm{2}}\left[(M_{\mathrm{N}}{x}_{\text{NG}}+M_{\mathrm{R}}{x}_{\text{NR}})\cos \left({\mathit{\nu }}_{\mathrm{y}}\right)\right.\\ \text{(43)}& {\displaystyle }& {\displaystyle }\left.+(M_{\mathrm{N}}{z}_{\text{NG}}+M_{\mathrm{R}}{z}_{\text{NR}})\sin \left({\mathit{\nu }}_{\mathrm{y}}q\right)\right]\text{(44)}& {\displaystyle }& {\displaystyle }+f_{\mathrm{a}}{\mathit{\nu }}_{\mathrm{y}}(x_{\text{NR}}\sin {\mathit{\theta }}_{\mathrm{t}}+z_{\text{NR}}\cos {\mathit{\theta }}_{\mathrm{t}})\text{(45)}& {\displaystyle }& {\displaystyle }+f_{\mathrm{a}}\cos \left({\mathit{\theta }}_{\mathrm{t}}+{\mathit{\nu }}_{\mathrm{y}}q\right).\end{array}$$

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ω²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_Nx_NG+M_Rx_NR) and (M_Nz_NG+M_Rz_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:

$$\begin{array}{}\text{(46)}& \begin{array}{rl}{K}_{\mathrm{q},\text{lin}}=& (K_{\mathrm{e}}+K_{\mathrm{g}})-{\mathit{\nu }}_{\mathrm{y}}^{\mathrm{2}}g(M_{\mathrm{N}}{z}_{\text{NG}}+M_{\mathrm{R}}{z}_{\text{NR}}\\ & -f_{\mathrm{a}}q\cos {\mathit{\theta }}_{\mathrm{t}})+{\mathit{\nu }}_{\mathrm{y}}{f}_{\mathrm{a}}\sin {\mathit{\theta }}_{\mathrm{t}}.\end{array}\end{array}$$

For the NREL 5 MW reference turbine (Jonkman et al., 2009), the different numerical values are g=9.807 m s⁻², θ_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², $J_{\oplus ,\mathrm{R}}=\mathrm{1.92}\times {\mathrm{10}}^{\mathrm{7}}$ kg m², $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², 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

$$\begin{array}{c}{\displaystyle }\begin{array}{rl}{\mathrm{\Phi }}_{\mathrm{1}}\left(z\right)=& \left(a_{\mathrm{2}}{\stackrel{\mathrm{‾}}{z}}^{\mathrm{2}}+a_{\mathrm{3}}{\stackrel{\mathrm{‾}}{z}}^{\mathrm{3}}+a_{\mathrm{4}}{\stackrel{\mathrm{‾}}{z}}^{\mathrm{4}}+a_{\mathrm{5}}{\stackrel{\mathrm{‾}}{z}}^{\mathrm{5}}+a_{\mathrm{6}}{\stackrel{\mathrm{‾}}{z}}^{\mathrm{6}}\right)/\\ & (a_{\mathrm{2}}+a_{\mathrm{3}}+a_{\mathrm{4}}+a_{\mathrm{5}}+a_{\mathrm{6}}),\end{array}\\ \text{(47)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{d{\mathrm{\Phi }}_{\mathrm{1}}}{\mathrm{d}z}}\left(z\right)=& {\displaystyle \frac{\mathrm{1}}{L_{\mathrm{T}}}}\left(\mathrm{2}{a}_{\mathrm{2}}\stackrel{\mathrm{‾}}{z}+\mathrm{3}{a}_{\mathrm{3}}{\stackrel{\mathrm{‾}}{z}}^{\mathrm{2}}+\mathrm{4}{a}_{\mathrm{4}}{\stackrel{\mathrm{‾}}{z}}^{\mathrm{3}}+\mathrm{5}{a}_{\mathrm{5}}{\stackrel{\mathrm{‾}}{z}}^{\mathrm{4}}\right.\\ & \left.+\mathrm{6}{a}_{\mathrm{6}}{\stackrel{\mathrm{‾}}{z}}^{\mathrm{5}}\right)/(a_{\mathrm{2}}+a_{\mathrm{3}}+a_{\mathrm{4}}+a_{\mathrm{5}}+a_{\mathrm{6}}),\end{array}{\displaystyle }\begin{array}{rl}{\displaystyle \frac{d^{\mathrm{2}}{\mathrm{\Phi }}_{\mathrm{1}}}{\mathrm{d}{z}^{\mathrm{2}}}}\left(z\right)=& {\displaystyle \frac{\mathrm{1}}{L_{\mathrm{T}}^{\mathrm{2}}}}\left(\mathrm{2}{a}_{\mathrm{2}}+\mathrm{6}{a}_{\mathrm{3}}\stackrel{\mathrm{‾}}{z}+\mathrm{12}{a}_{\mathrm{4}}{\stackrel{\mathrm{‾}}{z}}^{\mathrm{2}}+\mathrm{20}{a}_{\mathrm{5}}{\stackrel{\mathrm{‾}}{z}}^{\mathrm{3}}\right.\\ & \left.+\mathrm{30}{a}_{\mathrm{6}}{\stackrel{\mathrm{‾}}{z}}^{\mathrm{4}}\right)/(a_{\mathrm{2}}+a_{\mathrm{3}}+a_{\mathrm{4}}+a_{\mathrm{5}}+a_{\mathrm{6}}),\end{array}\end{array}$$

with $\stackrel{\mathrm{‾}}{z}=z/L$, a₂=0.7004, a₃=2.1963, $a_{\mathrm{4}}=-\mathrm{5.6202}$, a₅=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⁻¹. 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 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, Φ₁=Φ₁e_x and Φ₂=Φ₂e_y. The degrees of freedom are $\mathit{q}=(q_{\mathrm{1}},q_{\mathrm{2}},\mathit{\psi })$, where q₁ and q₂ 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 5 Properties of the NREL 5 MW turbine tower: mass per length (m), bending stiffness (EI), shape function displacement (Φ), slope ($d\mathrm{\Phi }/\mathrm{d}z$), and curvature $(d^{\mathrm{2}}\mathrm{\Phi }/\mathrm{d}{z}^{\mathrm{2}})$.

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 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 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², 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; 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.

Figure 6 Free 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 $\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 first-order small-angle approximation of θ_t and ϕ_y and a second-order approximation for ν_y.

Figure 7 Model 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=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.

Figure 8 Free-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.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 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 (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 low-order polynomial (usually linear or quadratic) approximations to eliminate high-order 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 “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 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:

$$\begin{array}{}\text{(A1)}& \begin{array}{rl}\left[\begin{array}{ccc}{\mathit{M}}_{\text{xx}}& {\mathit{M}}_{\mathrm{x}\mathit{\theta }}& {\mathit{M}}_{\text{xe}}\\ & {\mathit{M}}_{\mathit{\theta }\mathit{\theta }}& {\mathit{M}}_{\mathit{\theta }\mathrm{e}}\\ \text{sym.}& & {\mathit{M}}_{\text{ee}}\end{array}\right]\left[\begin{array}{c}{\mathit{a}}_{\mathrm{i}}\\ \dot{\mathit{\omega }}i\\ {\ddot{\mathit{q}}}_{\mathrm{e}}\end{array}\right]& +\left[\begin{array}{c}{\mathit{k}}_{\mathit{\omega },\mathrm{x}}\\ {\mathit{k}}_{\mathit{\omega },\mathit{\theta }}\\ {\mathit{k}}_{\mathit{\omega },\mathrm{e}}\end{array}\right]\\ & +\left[\begin{array}{c}\mathrm{0}\\ \mathrm{0}\\ {\mathit{k}}_{\mathrm{e}}\end{array}\right]=\left[\begin{array}{c}{\mathit{f}}_{\mathrm{x}}\\ {\mathit{f}}_{\mathit{\theta }}\\ {\mathit{f}}_{\mathrm{e}}\end{array}\right].\end{array}\end{array}$$

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

$$\begin{array}{}\text{(A2)}& {\displaystyle }{\mathit{M}}_{\text{xx}}& {\displaystyle }=\int {\mathit{I}}_{\mathrm{3}}\mathrm{d}m=M{\mathit{I}}_{\mathrm{3}}(\mathrm{3}\times \mathrm{3}),\text{(A3)}& {\displaystyle }{\mathit{M}}_{\mathrm{x}\mathit{\theta }}& {\displaystyle }=-\int {\stackrel{\mathrm{̃}}{\mathit{s}}}_{\mathrm{P}}\mathrm{d}m=-M{\stackrel{\mathrm{̃}}{\mathit{s}}}_{\text{CM}}(\mathrm{3}\times \mathrm{3}),\text{(A4)}& {\displaystyle }{\mathit{M}}_{\mathit{\theta }\mathit{\theta }}& {\displaystyle }=-\int {\stackrel{\mathrm{̃}}{\mathit{s}}}_{\mathrm{P}}{\stackrel{\mathrm{̃}}{\mathit{s}}}_{\mathrm{P}}\mathrm{d}m=\mathit{J}(\mathrm{3}\times \mathrm{3}),\text{(A5)}& {\displaystyle }{\mathit{M}}_{\mathit{\theta }\mathrm{e}}& {\displaystyle }=\int {\stackrel{\mathrm{̃}}{\mathit{s}}}_{\mathrm{P}}\mathbf{\Phi }\mathrm{d}m={\mathit{C}}_{\mathrm{r}}^{\mathrm{T}}(\mathrm{3}\times n_{\mathrm{e}}),\text{(A6)}& {\displaystyle }{\mathit{M}}_{\text{xe}}& {\displaystyle }=\int \mathbf{\Phi }\mathrm{d}m={\mathit{C}}_{\mathrm{t}}^{\mathrm{T}}(\mathrm{3}\times n_{\mathrm{e}}),\text{(A7)}& {\displaystyle }{\mathit{M}}_{\text{ee}}& {\displaystyle }=\int {\mathbf{\Phi }}^{\mathrm{T}}\mathbf{\Phi }\mathrm{d}m(n_{\mathrm{e}}\times n_{\mathrm{e}}).\end{array}$$

The integrals are volume integrals over the volume of the body (for beams, they reduce to line integrals). The notation $\left[\stackrel{\mathrm{̃}}{}\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 higher-order 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

$$\begin{array}{}\text{(A8)}& {\displaystyle }{\mathit{k}}_{\mathit{\omega },\mathrm{x}}=\mathrm{2}\stackrel{\mathrm{̃}}{\mathit{\omega }}{\mathit{C}}_{\mathrm{t}}^{\mathrm{T}}{\dot{\mathit{q}}}_{\mathrm{e}}+M\stackrel{\mathrm{̃}}{\mathit{\omega }}\stackrel{\mathrm{̃}}{\mathit{\omega }}{\mathit{s}}_{\text{CM}}(\mathrm{3}\times \mathrm{1}),\text{(A9)}& {\displaystyle }{\mathit{k}}_{\mathit{\omega },\mathit{\theta }}=\stackrel{\mathrm{̃}}{\mathit{\omega }}{\mathit{M}}_{\mathit{\theta }\mathit{\theta }}\mathit{\omega }+\left[\sum _{j=\mathrm{1}..n_{\mathrm{e}}}{\mathit{G}}_{\mathrm{r},j}{\dot{q}}_{\mathrm{e},j}\right]\mathit{\omega }(\mathrm{3}\times \mathrm{1}),\text{(A10)}& {\displaystyle }\begin{array}{rl}{\mathit{k}}_{\mathit{\omega },\mathrm{e}}=& {\left[{\mathit{\omega }}^{\mathrm{T}}{\mathit{O}}_{\mathrm{e},j}\mathit{\omega }\right]}_{j=\mathrm{1}..n_{\mathrm{e}}}\\ & +\left[\sum _{j=\mathrm{1}..n_{\mathrm{e}}}{\mathit{G}}_{\mathrm{e},j}{\dot{q}}_{\mathrm{e},j}\right]\mathit{\omega }(n_{\mathrm{e}}\times \mathrm{1}),\end{array}\end{array}$$

where

$$\begin{array}{}\text{(A11)}& {\displaystyle }{\mathit{G}}_{\mathrm{r},j}& {\displaystyle }=-\mathrm{2}\int {\stackrel{\mathrm{̃}}{\mathit{s}}}_{\mathrm{P}}{\stackrel{\mathrm{̃}}{\mathbf{\Phi }}}_j\mathrm{d}m(\mathrm{3}\times \mathrm{3}),\text{(A12)}& {\displaystyle }{\mathit{O}}_{\mathrm{e},j}& {\displaystyle }=\int {\stackrel{\mathrm{̃}}{\mathbf{\Phi }}}_j{\stackrel{\mathrm{̃}}{\mathit{s}}}_{\mathrm{P}}\mathrm{d}m=-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\mathit{G}}_{\mathrm{r},j}^{\mathrm{T}}(\mathrm{3}\times \mathrm{3}),\text{(A13)}& {\displaystyle }{\mathit{G}}_{\mathrm{e},j}& {\displaystyle }=-\mathrm{2}\int {\mathbf{\Phi }}^{\mathrm{T}}{\stackrel{\mathrm{̃}}{\mathbf{\Phi }}}_j\mathrm{d}m(n_{\mathrm{e}}\times \mathrm{3}).\end{array}$$

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

$$\begin{array}{}\text{(A14)}& {\displaystyle }\begin{array}{rl}{\mathit{O}}_{\mathrm{e}}=[& O_{\mathrm{e},j,\mathrm{11}},O_{\mathrm{e},j,\mathrm{22}},O_{\mathrm{e},j,\mathrm{33}},O_{\mathrm{e},j,\mathrm{12}}\\ & +O_{\mathrm{e},j,\mathrm{21}},O_{\mathrm{e},j,\mathrm{23}}+O_{\mathrm{e},j,\mathrm{32}},O_{\mathrm{e},j,\mathrm{13}}\\ & +O_{\mathrm{e},j,\mathrm{31}}{]}_{j=\mathrm{1}..n_{\mathrm{e}}}(n_{\mathrm{e}}\times \mathrm{6}),\end{array}\text{(A15)}& {\displaystyle }\mathbf{\Omega }=\left[{\mathit{\omega }}_{\mathrm{x}}^{\mathrm{2}},{\mathit{\omega }}_{\mathrm{y}}^{\mathrm{2}},{\mathit{\omega }}_{\mathrm{z}}^{\mathrm{2}},{\mathit{\omega }}_{\mathrm{x}}{\mathit{\omega }}_{\mathrm{y}},{\mathit{\omega }}_{\mathrm{y}}{\mathit{\omega }}_{\mathrm{z}},{\mathit{\omega }}_{\mathrm{x}}{\mathit{\omega }}_{\mathrm{z}}\right](\mathrm{6}\times \mathrm{1}).\end{array}$$

The body elastic forces are given by

$$\begin{array}{}\text{(A16)}& {\displaystyle }{\mathit{k}}_{\mathrm{e}}={\mathit{k}}_{\mathit{\sigma }}+{\mathit{K}}_{\mathrm{e}}{\mathit{q}}_{\mathrm{e}}+{\mathit{D}}_{\mathrm{e}}{\dot{\mathit{q}}}_{\mathrm{e}},\end{array}$$

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:

$$\begin{array}{}\text{(A17)}& {\displaystyle }{\mathit{f}}_{\mathrm{x}}& {\displaystyle }=\int \mathit{p}\mathrm{d}V+{\mathit{M}}_{\text{xx}}\mathit{g}(\mathrm{3}\times \mathrm{1}),\text{(A18)}& {\displaystyle }{\mathit{f}}_{\mathit{\theta }}& {\displaystyle }=\int {\mathit{s}}_{\mathrm{P}}\times \mathit{p}\mathrm{d}V+{\mathit{M}}_{\mathit{\theta }x}\mathit{g}(\mathrm{3}\times \mathrm{1}),\text{(A19)}& {\displaystyle }{\mathit{f}}_{\mathrm{e}}& {\displaystyle }=\int {\mathbf{\Phi }}^{\mathrm{T}}\mathit{p}\mathrm{d}V+{\mathit{M}}_{\text{ex}}\mathit{g}(n_{\mathrm{e}}\times \mathrm{1}).\end{array}$$

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, 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)

$$\begin{array}{}\text{(B1)}& {\displaystyle }{M}_{\mathrm{e},ij}=\underset{\mathrm{0}}{\overset{L}{\int }}m\left(z\right){\mathbf{\Phi }}_i\left(z\right)\cdot {\mathbf{\Phi }}_j\left(z\right)\mathrm{d}z,i,j=\mathrm{1},\mathrm{\dots }{n}_{\mathrm{e}}.\end{array}$$

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

$$\begin{array}{}\text{(B2)}& \begin{array}{rl}{K}_{\mathrm{e},ij}=\underset{\mathrm{0}}{\overset{L}{\int }}[& {\text{EI}}_{\mathrm{y}}{\displaystyle \frac{d^{\mathrm{2}}{\mathrm{\Phi }}_{i,\mathrm{x}}}{\mathrm{d}{z}^{\mathrm{2}}}}{\displaystyle \frac{d^{\mathrm{2}}{\mathrm{\Phi }}_{j,\mathrm{x}}}{\mathrm{d}{z}^{\mathrm{2}}}}\\ & +{\text{EI}}_{\mathrm{x}}{\displaystyle \frac{d^{\mathrm{2}}{\mathrm{\Phi }}_{i,y}}{\mathrm{d}{z}^{\mathrm{2}}}}{\displaystyle \frac{d^{\mathrm{2}}{\mathrm{\Phi }}_{j,y}}{\mathrm{d}{z}^{\mathrm{2}}}}]\mathrm{d}z,i,j=\mathrm{1},\mathrm{\dots }{n}_{\mathrm{e}}.\end{array}\end{array}$$

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)

$$\begin{array}{}\text{(B3)}& {\displaystyle }{f}_{\mathrm{e},\mathrm{i}}=\underset{\mathrm{0}}{\overset{L}{\int }}{\mathbf{\Phi }}_i\cdot \mathit{p}\mathrm{d}z.\end{array}$$

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

$$\begin{array}{}\text{(B4)}& {\displaystyle }{\mathit{M}}_{\mathrm{e}}{\ddot{\mathit{q}}}_{\mathrm{e}}+{\mathit{D}}_{\mathrm{e}}{\dot{\mathit{q}}}_{\mathrm{e}}+{\mathit{K}}_{\mathrm{e}}{\mathit{q}}_{\mathrm{e}}={\mathit{f}}_{\mathrm{e}},\end{array}$$

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}\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

$$\begin{array}{}\text{(B5)}& \begin{array}{rl}{K}_{\mathrm{N},ij}& =-\underset{\mathrm{0}}{\overset{L}{\int }}{\mathbf{\Phi }}_i\cdot {\displaystyle \frac{\mathrm{d}}{\mathrm{d}z}}\left[N\left(z\right){\displaystyle \frac{\mathrm{d}{\mathbf{\Phi }}_{\mathit{j}}}{\mathrm{d}z}}\right]\mathrm{d}z\\ & =\underset{\mathrm{0}}{\overset{L}{\int }}N\left(z\right){\displaystyle \frac{\mathrm{d}{\mathbf{\Phi }}_i}{\mathrm{d}z}}\cdot {\displaystyle \frac{\mathrm{d}{\mathbf{\Phi }}_j}{\mathrm{d}z}}-{\left[N\left(z\right){\mathbf{\Phi }}_{\mathbf{i}}\cdot {\displaystyle \frac{\mathrm{d}{\mathbf{\Phi }}_{\mathit{j}}}{\mathrm{d}z}}\right]}_{\mathrm{0}}^L,\end{array}\end{array}$$

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

$$\begin{array}{}\text{(C1)}& {\displaystyle }{\mathit{k}}_{\mathit{\sigma }}={\mathit{K}}_{\mathrm{g}}{\mathit{q}}_{\mathrm{e}},\end{array}$$

where K_g is the geometric stiffness matrix of shape n_e×n_e. 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 ($\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:

$$\begin{array}{}\text{(C2)}& \begin{array}{rl}{\mathit{K}}_{\mathrm{g}}=& \sum _{\mathit{\alpha }=\mathrm{1}}^{\mathrm{3}}\left[\left(a_{\mathit{\alpha }}-g_{\mathit{\alpha }}\right){\mathit{K}}_{\text{gt},\mathit{\alpha }}+{\dot{\mathit{\omega }}}_{\mathit{\alpha }}{\mathit{K}}_{\text{gr},\mathit{\alpha }}\right]\\ & +\sum _{\mathit{\alpha }=\mathrm{1}}^{\mathrm{3}}\sum _{\mathit{\beta }=\mathrm{1}}^{\mathrm{3}}{\mathit{\omega }}_{\mathit{\alpha }}{\mathit{\omega }}_{\mathit{\beta }}{\mathit{K}}_{\mathrm{g}\mathit{\omega },\mathit{\alpha }\mathit{\beta }}\\ & +\sum _{\mathit{\alpha }=\mathrm{1}}^{\mathrm{3}}\left[p_{\mathit{\alpha }}{\mathit{K}}_{\text{gp},\mathit{\alpha }}+\sum _{\mathrm{k}}\left(F_{\mathit{\alpha }}^k{\mathit{K}}_{\text{gF},\mathit{\alpha }}^k+{\mathit{\tau }}_{\mathit{\alpha }}^k{\mathit{K}}_{\mathrm{g}\mathit{\tau },\mathit{\alpha }}^k\right)\right],\end{array}\end{array}$$

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 steady-state 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² m⁻² .

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:

$$\begin{array}{}\text{(C3)}& {\displaystyle }{k}_{\mathit{\sigma },N,\mathrm{i}}=\underset{\mathrm{0}}{\overset{L}{\int }}N\left(z\right){\mathbf{\Phi }}_i\cdot \left[\sum _j{\displaystyle \frac{\mathrm{d}{\mathbf{\Phi }}_j}{\mathrm{d}z}}{q}_j\right]\mathrm{d}z.\end{array}$$

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:

$$\begin{array}{}\text{(C4)}& \begin{array}{rl}& K_{\text{gt},z,ij}=\underset{\mathrm{0}}{\overset{L}{\int }}N\left(z\right){\displaystyle \frac{\mathrm{d}{\mathbf{\Phi }}_{\mathbf{i}}}{\mathrm{d}z}}\cdot {\displaystyle \frac{\mathrm{d}{\mathbf{\Phi }}_{\mathit{j}}}{\mathrm{d}z}}\mathrm{d}z,\\ & N\left(z\right)=\underset{\mathrm{z}}{\overset{L}{\int }}m\left(z\right)\mathrm{d}z.\end{array}\end{array}$$

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

$$\begin{array}{}\text{(C5)}& \begin{array}{rl}& K_{\text{gF},z,ij}^k=\underset{\mathrm{0}}{\overset{L}{\int }}N\left(z\right){\displaystyle \frac{\mathrm{d}{\mathbf{\Phi }}_{\mathbf{i}}}{\mathrm{d}z}}\cdot {\displaystyle \frac{\mathrm{d}{\mathbf{\Phi }}_{\mathit{j}}}{\mathrm{d}z}}\mathrm{d}z,\\ & N\left(z\right)=\mathrm{1}\text{ if }z<z_{\mathrm{k}},\mathrm{0}\text{ otherwise}.\end{array}\end{array}$$

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:

$$\begin{array}{}\text{(C6)}& \begin{array}{rl}& K_{\mathrm{g}\mathit{\omega },\mathit{\alpha }\mathit{\beta },ij}=\underset{\mathrm{0}}{\overset{L}{\int }}-{\mathit{e}}_{\mathrm{z}}\cdot \left({\stackrel{\mathrm{̃}}{\mathit{e}}}_{\mathit{\alpha }}{\stackrel{\mathrm{̃}}{\mathit{e}}}_{\mathit{\beta }}\mathit{N}\left(z\right)\right){\displaystyle \frac{\mathrm{d}{\mathbf{\Phi }}_{\mathbf{i}}}{\mathrm{d}z}}\cdot {\displaystyle \frac{\mathrm{d}{\mathbf{\Phi }}_{\mathit{j}}}{\mathrm{d}z}}\mathrm{d}z,\\ & \mathit{N}\left(z\right)=\underset{\mathrm{z}}{\overset{L}{\int }}m\left(z\right){\mathit{s}}_{P_{\mathrm{0}}}\mathrm{d}z.\end{array}\end{array}$$

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

$$\begin{array}{}\text{(C7)}& \begin{array}{rl}& K_{\text{gr},\mathit{\alpha },ij}=\underset{\mathrm{0}}{\overset{L}{\int }}-{\mathit{e}}_{\mathrm{z}}\cdot \left({\stackrel{\mathrm{̃}}{\mathit{e}}}_{\mathit{\alpha }}\mathit{N}\left(z\right)\right){\displaystyle \frac{\mathrm{d}{\mathbf{\Phi }}_{\mathbf{i}}}{\mathrm{d}z}}\cdot {\displaystyle \frac{\mathrm{d}{\mathbf{\Phi }}_{\mathit{j}}}{\mathrm{d}z}}\mathrm{d}z,\\ & \mathit{N}\left(z\right)=\underset{\mathrm{z}}{\overset{L}{\int }}m\left(z\right){\mathit{s}}_{P_{\mathrm{0}}}\mathrm{d}z.\end{array}\end{array}$$

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

$$\begin{array}{}\text{(C8)}& \begin{array}{rl}& k_{\mathit{\sigma },\text{Cor},\mathrm{i}}=\underset{\mathrm{0}}{\overset{L}{\int }}N\left(z\right){\mathbf{\Phi }}_i\cdot \left[\sum _j{\displaystyle \frac{\mathrm{d}{\mathbf{\Phi }}_j}{\mathrm{d}z}}{q}_j\right]\mathrm{d}z,\\ & N\left(z\right)=\mathrm{2}\underset{\mathrm{z}}{\overset{L}{\int }}m\left(z\right)\sum _k\left[{\mathit{\omega }}_{\mathrm{x}}{\mathrm{\Phi }}_{\mathrm{k},y}\left(z\right)-{\mathit{\omega }}_{\mathrm{y}}{\mathrm{\Phi }}_{\mathrm{k},x}\left(z\right)\right]{\dot{q}}_{\mathrm{k}}\mathrm{d}z.\end{array}\end{array}$$

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 ${\mathit{k}}_{\mathit{\sigma },\text{Cor}}({\mathit{q}}_{\mathrm{e}},{\dot{\mathit{q}}}_{\mathrm{e}})$ is added to the right-hand side of Eq. (C1).

C3 Integration into the equations of motion

The term k_σ=K_gq_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 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., Lemmer, 2018, and Branlard, 2019), the Jacobians are sometimes stacked into a matrix form:

$$\begin{array}{}\text{(D1)}& {\displaystyle }\mathit{J}=\left[\begin{array}{c}{\mathit{J}}_{\mathrm{v}}\\ {\mathit{J}}_{\mathit{\omega }}\\ {\mathit{J}}_{\mathrm{e}}\end{array}\right].\end{array}$$

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 (Wallrapp, 1994; Schwertassek and Wallrapp, 1999; Lemmer, 2018). In such a formulation, the linear rotation field is obtained as $\mathit{I}+\stackrel{\mathrm{̃}}{\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 time-varying 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

$$\begin{array}{}\text{(D2)}& \mathit{T}\left({\mathit{q}}_{\mathrm{e}}\right)={\mathit{T}}^{\mathrm{0}}+\sum _{j=\mathrm{1}..n_{\mathrm{e}}}{\mathit{T}}_j^{\mathrm{1}}{q}_{\mathrm{e},j}.\end{array}$$

If T is an array, T⁰ 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

$$\begin{array}{}\text{(D3)}& {\mathit{M}}_{\mathrm{x}\mathit{\theta }}=-\int {\stackrel{\mathrm{̃}}{\mathit{s}}}_{\mathrm{P}}\mathrm{d}m={\mathit{M}}_{\mathrm{x}\mathit{\theta }}^{\mathrm{0}}+\sum _{j=\mathrm{1}..n_{\mathrm{e}}}{\mathit{M}}_{\mathrm{x}\mathit{\theta },j}^{\mathrm{1}}{q}_{\mathrm{e},j},\end{array}$$

with

$$\begin{array}{}\text{(D4)}& {\mathit{M}}_{\mathrm{x}\mathit{\theta }}^{\mathrm{0}}=-\int {\stackrel{\mathrm{̃}}{\mathit{s}}}_{P_{\mathrm{0}}}\mathrm{d}m,{\mathit{M}}_{\mathrm{x}\mathit{\theta },j}^{\mathrm{1}}=-\int {\stackrel{\mathrm{̃}}{\mathbf{\Phi }}}_j\mathrm{d}m.\end{array}$$

The zeroth- and first-order 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 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 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 first-order terms of Eq. (D2). For instance, if we write M_ex (given in Eq. A6) using a first-order 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 first-order 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 Φ¹, and K_gτ into Ψ¹ 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 ${\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 Φ⁰ and Φ¹ 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 Φ¹. The advantage of this method is that the Φ¹ terms effectively account for the geometric stiffness. In practice, it is equivalent, and as convenient, to neglect the Φ¹ 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:

$$\begin{array}{}\text{(D5)}& f_r+f_r^{\ast }=\mathrm{0}={\mathit{J}}_{v_{\mathrm{P}},r}\cdot {\mathit{f}}_{\mathrm{P}}+{\mathit{J}}_{{\mathit{\omega }}_{\mathrm{P}},r}\cdot {\mathit{\tau }}_{\mathrm{P}}+h_{\mathrm{r}},\end{array}$$

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”):

$$\begin{array}{}\text{(D6)}& {\mathit{f}}_{\mathrm{P}}=\sum _{j=\mathrm{1}}^{n_{\mathrm{q}}}{\mathit{f}}_{\mathrm{P},j}{\ddot{q}}_j+{\mathit{f}}_{\mathrm{P},t},{\mathit{\tau }}_{\mathrm{P}}=\sum _{j=\mathrm{1}}^{n_{\mathrm{q}}}{\mathit{\tau }}_{\mathrm{P},j}{\ddot{q}}_j+{\mathit{\tau }}_{\mathrm{P},t}.\end{array}$$

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 right-hand side of the equation of motion (Eq. 23) are obtained as

$$\begin{array}{}\text{(D7)}& \begin{array}{rl}& M_{\text{rj}}=-{\mathit{J}}_{v_{\mathrm{P}},r}\cdot {\mathit{f}}_{\mathrm{P},j}-{\mathit{J}}_{{\mathit{\omega }}_{\mathrm{P}},r}\cdot {\mathit{\tau }}_{\mathrm{P},j},\\ & F_{\mathrm{r}}={\mathit{J}}_{v_{\mathrm{P}},r}\cdot {\mathit{f}}_{\mathrm{P},t}+{\mathit{J}}_{{\mathit{\omega }}_{\mathrm{P}},r}\cdot {\mathit{\tau }}_{\mathrm{P},t}+h_{\mathrm{r}}.\end{array}\end{array}$$

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, $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:

$$\begin{array}{}\text{(E1)}& {\mathit{v}}_{\mathrm{O},\mathrm{T}}=\mathbf{0},{\mathit{\omega }}_{\mathrm{T}}=\mathbf{0},{\mathit{a}}_{\mathrm{O},\mathrm{T}}=\mathbf{0}.\end{array}$$

All Jacobians are zero except ${\mathit{J}}_{\mathrm{e},\mathrm{1}T}=\mathrm{1}$. The inertial force, torque, and elastic force are

$$\begin{array}{}\text{(E2)}& \begin{array}{rl}& {\mathit{f}}_{\mathrm{T}}^{\ast }=C_{tTx}\ddot{q}{\widehat{\mathit{t}}}_{\mathrm{x}}+M_{\mathrm{T}}g{\widehat{\mathit{t}}}_{\mathrm{z}},{\mathit{\tau }}_{\mathrm{T}}^{\ast }=C_{rTy}\ddot{q}{\widehat{\mathit{t}}}_{\mathrm{y}},\\ & {\mathit{E}}_{\mathrm{T}}^{\ast }=f_{\mathrm{e}}+D_{\mathrm{e}}\dot{q}+(K_{\mathrm{e}}+K_{\mathrm{q}})q+M_{\mathrm{e}}\ddot{q}.\end{array}\end{array}$$

The nacelle kinematics (at its center of mass) are

$$\begin{array}{}\text{(E3)}& {\displaystyle }{\mathit{v}}_{\mathrm{G},N}=\dot{q}{\widehat{\mathit{t}}}_{\mathrm{x}}+{\mathit{\nu }}_{\mathrm{y}}{z}_{\text{NG}}\dot{q}{\widehat{\mathit{n}}}_{\mathrm{x}}-{\mathit{\nu }}_{\mathrm{y}}{x}_{\text{NG}}\dot{q}{\widehat{\mathit{n}}}_{\mathrm{z}},{\mathit{\omega }}_{\mathrm{N}}={\mathit{\nu }}_{\mathrm{y}}\dot{q}{\widehat{\mathit{t}}}_{\mathrm{y}},\text{(E4)}& {\displaystyle }\begin{array}{rl}{\mathit{a}}_{\mathrm{G},N}=& \ddot{q}{\widehat{\mathit{t}}}_{\mathrm{x}}+(-{\mathit{\nu }}_{\mathrm{y}}^{\mathrm{2}}{x}_{\text{NG}}{\dot{q}}^{\mathrm{2}}+{\mathit{\nu }}_{\mathrm{y}}{z}_{\text{NG}}\ddot{q}){\widehat{\mathit{n}}}_{\mathrm{x}}\\ & +(-{\mathit{\nu }}_{\mathrm{y}}^{\mathrm{2}}{z}_{\text{NG}}{\dot{q}}^{\mathrm{2}}-{\mathit{\nu }}_{\mathrm{y}}{x}_{\text{NG}}\ddot{q}){\widehat{\mathit{n}}}_{\mathrm{z}}.\end{array}\end{array}$$

The Jacobians with respect to q are

$$\begin{array}{}\text{(E5)}& {\mathit{J}}_{\mathrm{v},\mathrm{1}N}={\widehat{\mathit{t}}}_{\mathrm{x}}+{\mathit{\nu }}_{\mathrm{y}}{z}_{\text{NG}}{\widehat{\mathit{n}}}_{\mathrm{x}}-{\mathit{\nu }}_{\mathrm{y}}{x}_{\text{NG}}{\widehat{\mathit{n}}}_{\mathrm{z}},{\mathit{J}}_{\mathit{\omega },\mathrm{1}N}={\mathit{\nu }}_{\mathrm{y}}{\widehat{\mathit{t}}}_{\mathrm{y}}.\end{array}$$

The inertial force and torque on the nacelle are

$$\begin{array}{}\text{(E6)}& \begin{array}{rl}& {\mathit{f}}_{\mathrm{N}}^{\ast }=M_{\mathrm{N}}\ddot{q}{\widehat{\mathit{t}}}_{\mathrm{x}}+M_{\mathrm{N}}\left(-{\mathit{\nu }}_{\mathrm{y}}^{\mathrm{2}}{x}_{\text{NG}}{\dot{q}}^{\mathrm{2}}+{\mathit{\nu }}_{\mathrm{y}}{z}_{\text{NG}}\ddot{q}\right){\widehat{\mathit{n}}}_{\mathrm{x}}\\ & +M_{\mathrm{N}}\left(-{\mathit{\nu }}_{\mathrm{y}}^{\mathrm{2}}{z}_{\text{NG}}{\dot{q}}^{\mathrm{2}}-{\mathit{\nu }}_{\mathrm{y}}{x}_{\text{NG}}\ddot{q}\right){\widehat{\mathit{n}}}_{\mathrm{z}},\\ & {\mathit{\tau }}_{\mathrm{N}}^{\ast }=J_{\text{y,N}}{\mathit{\nu }}_{\mathrm{y}}\ddot{q}{\widehat{\mathit{n}}}_{\mathrm{y}}.\end{array}\end{array}$$

The kinematics of the rotor are

$$\begin{array}{}\text{(E7)}& {\displaystyle }\begin{array}{rl}& {\mathit{v}}_{\mathrm{G},R}=\dot{q}{\widehat{\mathit{t}}}_{\mathrm{x}}+{\mathit{\nu }}_{\mathrm{y}}{z}_{\text{NR}}\dot{q}{\widehat{\mathit{n}}}_{\mathrm{x}}-{\mathit{\nu }}_{\mathrm{y}}{x}_{\text{NR}}\dot{q}{\widehat{\mathit{n}}}_{\mathrm{z}},\\ & {\mathit{\omega }}_{\mathrm{R}}=\dot{\mathit{\psi }}{\widehat{\mathit{e}}}_{{\mathbf{r}}_{\mathbf{x}}}+{\mathit{\nu }}_{\mathrm{y}}\dot{q}{\widehat{\mathit{t}}}_{\mathrm{y}},\end{array}\text{(E8)}& {\displaystyle }\begin{array}{rl}{\mathit{a}}_{\mathrm{G},R}=& \ddot{q}{\widehat{\mathit{t}}}_{\mathrm{x}}+(-{\mathit{\nu }}_{\mathrm{y}}^{\mathrm{2}}{x}_{\text{NR}}{\dot{q}}^{\mathrm{2}}+{\mathit{\nu }}_{\mathrm{y}}{z}_{\text{NR}}\ddot{q}){\widehat{\mathit{n}}}_{\mathrm{x}}\\ & +(-{\mathit{\nu }}_{\mathrm{y}}^{\mathrm{2}}{z}_{\text{NR}}{\dot{q}}^{\mathrm{2}}-{\mathit{\nu }}_{\mathrm{y}}{x}_{\text{NR}}\ddot{q}){\widehat{\mathit{n}}}_{\mathrm{z}}.\end{array}\end{array}$$

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

$$\begin{array}{rl}& {\mathit{J}}_{\mathrm{v},\mathrm{1}R}={\widehat{\mathit{t}}}_{\mathrm{x}}+{\mathit{\nu }}_{\mathrm{y}}{z}_{\text{NR}}{\widehat{\mathit{n}}}_{\mathrm{x}}-{\mathit{\nu }}_{\mathrm{y}}{x}_{\text{NR}}{\widehat{\mathit{n}}}_{\mathrm{z}},{\mathit{J}}_{\mathit{\omega },\mathrm{1}R}={\mathit{\nu }}_{\mathrm{y}}{\widehat{\mathit{t}}}_{\mathrm{y}},\\ & {\mathit{J}}_{\mathit{\omega },\mathrm{2}R}={\widehat{\mathit{r}}}_{\mathrm{x}}.\end{array}$$

The inertial force and torque on the rotor are

$$\begin{array}{ll}{\displaystyle }{\mathit{f}}_{\mathrm{R}}^{\ast }=& {\displaystyle }{M}_{\mathrm{R}}\ddot{q}{\widehat{\mathit{t}}}_{\mathrm{x}}+M_{\mathrm{R}}\left(-{\mathit{\nu }}_{\mathrm{y}}^{\mathrm{2}}{x}_{\text{NR}}{\dot{q}}^{\mathrm{2}}+{\mathit{\nu }}_{\mathrm{y}}{z}_{\text{NR}}\ddot{q}\right){\widehat{\mathit{n}}}_{\mathrm{x}}\\ \text{(E9)}& {\displaystyle }& {\displaystyle }+M_{\mathrm{R}}\left(-{\mathit{\nu }}_{\mathrm{y}}^{\mathrm{2}}{z}_{\text{NR}}{\dot{q}}^{\mathrm{2}}-{\mathit{\nu }}_{\mathrm{y}}{x}_{\text{NR}}\ddot{q}\right){\widehat{\mathit{n}}}_{\mathrm{z}},\text{(E10)}& {\displaystyle }{\mathit{\tau }}_{\mathrm{R}}^{\ast }=& {\displaystyle }{J}_{\text{x,R}}\ddot{\mathit{\psi }}{\widehat{\mathit{r}}}_{\mathrm{x}}{\displaystyle }& {\displaystyle }+\left(J_{\oplus ,\mathrm{R}}{\mathit{\nu }}_{\mathrm{y}}\sin \right(\mathit{\psi })\dot{\mathit{\psi }}\dot{q}+J_{\oplus ,\mathrm{R}}\left(-{\mathit{\nu }}_{\mathrm{y}}\sin \left(\mathit{\psi }\right)\dot{\mathit{\psi }}\dot{q}\right.\\ \text{(E11)}& {\displaystyle }& {\displaystyle }\left.+{\mathit{\nu }}_{\mathrm{y}}\cos \left(\mathit{\psi }\right)\ddot{q}\right)-J_{\text{x,R}}{\mathit{\nu }}_{\mathrm{y}}\sin \left(\mathit{\psi }\right)\dot{\mathit{\psi }}\dot{q}){\widehat{\mathit{r}}}_{\mathrm{y}}{\displaystyle }& {\displaystyle }+\left(J_{\oplus ,\mathrm{R}}{\mathit{\nu }}_{\mathrm{y}}\cos \right(\mathit{\psi })\dot{\mathit{\psi }}\dot{q}+J_{\oplus ,\mathrm{R}}\left(-{\mathit{\nu }}_{\mathrm{y}}\sin \left(\mathit{\psi }\right)\ddot{q}\right.\\ \text{(E12)}& {\displaystyle }& {\displaystyle }\left.-{\mathit{\nu }}_{\mathrm{y}}\cos \left(\mathit{\psi }\right)\dot{\mathit{\psi }}\dot{q}\right)-J_{\text{x,R}}{\mathit{\nu }}_{\mathrm{y}}\cos \left(\mathit{\psi }\right)\dot{\mathit{\psi }}\dot{q}){\widehat{\mathit{r}}}_{\mathrm{z}}.\end{array}$$

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 $\mathit{q}=(q_{\mathrm{1}},q_{\mathrm{2}},\mathit{\psi })$, are given in this section. The elements of the mass matrix are

$$\begin{array}{ll}\text{(E13)}& {\displaystyle }{M}_{\mathrm{11}}=& {\displaystyle }\left[M_{\mathrm{e}\mathrm{11}}+M_{\mathrm{N}}+M_{\mathrm{R}}\right]{\displaystyle }& {\displaystyle }+\left[J_{\text{y,N}}+J_{\oplus ,\mathrm{R}}+M_{\mathrm{N}}\left(x_{\text{NG}}^{\mathrm{2}}-\mathrm{2}{x}_{\text{NG}}{q}_{\mathrm{1}}+z_{\text{NG}}^{\mathrm{2}}\right)\right.\\ \text{(E14)}& {\displaystyle }& {\displaystyle }\left.+M_{\mathrm{R}}\left(x_{\text{NR}}^{\mathrm{2}}-\mathrm{2}{x}_{\text{NR}}{q}_{\mathrm{1}}+z_{\text{NR}}^{\mathrm{2}}\right)\right]{\mathit{\nu }}_{\mathrm{y}}^{\mathrm{2}}\text{(E15)}& {\displaystyle }& {\displaystyle }+\mathrm{2}\left[M_{\mathrm{N}}{z}_{\text{NG}}+M_{\mathrm{R}}{z}_{\text{NR}}\right]{\mathit{\nu }}_{\mathrm{y}},\text{(E16)}& {\displaystyle }{M}_{\mathrm{13}}=& {\displaystyle }{J}_{\text{x,R}}{\mathit{\theta }}_{\mathrm{t}}{\mathit{\nu }}_{\mathrm{x}}{\mathit{\nu }}_{\mathrm{y}}{q}_{\mathrm{2}},\text{(E17)}& {\displaystyle }{M}_{\mathrm{22}}=& {\displaystyle }\left[M_{\mathrm{e}\mathrm{22}}+M_{\mathrm{N}}+M_{\mathrm{R}}\right]\text{(E18)}& {\displaystyle }& {\displaystyle }+\left[J_{\text{x,N}}+J_{\text{x,R}}+M_{\mathrm{N}}{z}_{\text{NG}}^{\mathrm{2}}+M_{\mathrm{R}}{z}_{\text{NR}}^{\mathrm{2}}\right]{\mathit{\nu }}_{\mathrm{x}}^{\mathrm{2}}\text{(E19)}& {\displaystyle }& {\displaystyle }-\mathrm{2}\left[M_{\mathrm{N}}{z}_{\text{NG}}+M_{\mathrm{R}}{z}_{\text{NR}}\right]{\mathit{\nu }}_{\mathrm{x}},\text{(E20)}& {\displaystyle }{M}_{\mathrm{23}}=& {\displaystyle }{J}_{\text{x,R}}{\mathit{\nu }}_{\mathrm{x}},\text{(E21)}& {\displaystyle }{M}_{\mathrm{33}}=& {\displaystyle }{J}_{\text{x,R}}.\end{array}$$

The elements of the forcing vector are

$$\begin{array}{ll}{\displaystyle }{f}_{\mathrm{1}}=& {\displaystyle }{f}_{\mathrm{e}\mathrm{1}}-K_{\mathrm{e}\mathrm{11}}{q}_{\mathrm{1}}-D_{\mathrm{e}\mathrm{11}}{\dot{q}}_{\mathrm{1}}-J_{\text{x,R}}{\mathit{\theta }}_{\mathrm{t}}{\mathit{\nu }}_{\mathrm{x}}{\mathit{\nu }}_{\mathrm{y}}\dot{\mathit{\psi }}{\dot{q}}_{\mathrm{2}}\\ \text{(E22)}& {\displaystyle }& {\displaystyle }+\left[M_{\mathrm{N}}{x}_{\text{NG}}+M_{\mathrm{R}}{x}_{\text{NR}}\right]{\mathit{\nu }}_{\mathrm{y}}^{\mathrm{2}}{\dot{q}}_{\mathrm{1}}^{\mathrm{2}}{\displaystyle }& {\displaystyle }+g\left[M_{\mathrm{N}}\left({\mathit{\nu }}_{\mathrm{y}}^{\mathrm{2}}{z}_{\text{NG}}{q}_{\mathrm{1}}+{\mathit{\nu }}_{\mathrm{y}}{x}_{\text{NG}}\right)\right.\\ {\displaystyle }& {\displaystyle }\left.+M_{\mathrm{R}}\left({\mathit{\nu }}_{\mathrm{y}}^{\mathrm{2}}{z}_{\text{NR}}{q}_{\mathrm{1}}+{\mathit{\nu }}_{\mathrm{y}}{x}_{\text{NR}}\right)\right]\\ \text{(E23)}& {\displaystyle }& {\displaystyle }+f_{\mathrm{a}}\left[{\mathit{\theta }}_{\mathrm{t}}{\mathit{\nu }}_{\mathrm{y}}{x}_{\text{NR}}-{\mathit{\theta }}_{\mathrm{t}}{\mathit{\nu }}_{\mathrm{y}}{q}_{\mathrm{1}}+{\mathit{\nu }}_{\mathrm{y}}{z}_{\text{NR}}+\mathrm{1}\right],\text{(E24)}& {\displaystyle }{f}_{\mathrm{2}}=& {\displaystyle }{f}_{\mathrm{e}\mathrm{2}}-K_{\mathrm{e}\mathrm{22}}{q}_{\mathrm{2}}-D_{\mathrm{e}\mathrm{22}}{\dot{q}}_{\mathrm{2}}+J_{\text{x,R}}{\mathit{\theta }}_{\mathrm{t}}{\mathit{\nu }}_{\mathrm{x}}{\mathit{\nu }}_{\mathrm{y}}\dot{\mathit{\psi }}{\dot{q}}_{\mathrm{1}}\text{(E25)}& {\displaystyle }& {\displaystyle }+g\left[M_{\mathrm{N}}{z}_{\text{NG}}+M_{\mathrm{R}}{z}_{\text{NR}}\right]{\mathit{\nu }}_{\mathrm{x}}^{\mathrm{2}}{q}_{\mathrm{2}}+f_{\mathrm{a}}{\mathit{\theta }}_{\mathrm{t}}{\mathit{\nu }}_{\mathrm{x}}{q}_{\mathrm{2}},\text{(E26)}& {\displaystyle }{f}_{\mathrm{3}}=& {\displaystyle }-J_{\text{x,R}}{\mathit{\theta }}_{\mathrm{t}}{\mathit{\nu }}_{\mathrm{x}}{\mathit{\nu }}_{\mathrm{y}}{\dot{q}}_{\mathrm{1}}{\dot{q}}_{\mathrm{2}}+{\mathit{\tau }}_{\mathrm{a}}.\end{array}$$

E3 Three-degrees-of-freedom 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

$$\begin{array}{ll}\text{(E27)}& {\displaystyle }{M}_{\mathrm{11}}=& {\displaystyle }{M}_{\mathrm{F}}+M_{\mathrm{T}}+M_{\mathrm{N}},{\displaystyle }{M}_{\mathrm{12}}=& {\displaystyle }{M}_{\mathrm{F}}{z}_{\text{FG}}-M_{dTz}+M_{\mathrm{N}}\left[L_{\mathrm{T}}+z_{\text{NG}}\right.\\ {\displaystyle }& {\displaystyle }-{\mathit{\nu }}_{\mathrm{y}}{x}_{\text{NG}}{q}_{\mathrm{T}}\\ \text{(E28)}& {\displaystyle }& {\displaystyle }\left.-{\mathit{\varphi }}_{\mathrm{y}}(x_{\text{NG}}+q_{\mathrm{T}}+{\mathit{\nu }}_{\mathrm{y}}{z}_{\text{NG}}{q}_{\mathrm{T}})\right],{\displaystyle }{M}_{\mathrm{13}}=& {\displaystyle }{C}_{tT\mathrm{1}x}+M_{\mathrm{N}}\left[\mathrm{1}+{\mathit{\nu }}_{\mathrm{y}}{z}_{\text{NG}}-{\mathit{\nu }}_{\mathrm{y}}^{\mathrm{2}}{x}_{\text{NG}}{q}_{\mathrm{T}}\right.\\ \text{(E29)}& {\displaystyle }& {\displaystyle }\left.-{\mathit{\varphi }}_{\mathrm{y}}({\mathit{\nu }}_{\mathrm{y}}^{\mathrm{2}}{z}_{\text{NG}}{q}_{\mathrm{T}}+{\mathit{\nu }}_{\mathrm{y}}{x}_{\text{NG}})\right],{\displaystyle }{M}_{\mathrm{22}}=& {\displaystyle }{J}_{y,F}+M_{\mathrm{F}}{z}_{\text{FG}}^{\mathrm{2}}+J_{\mathrm{T},y}+J_{\text{y,N}}\\ {\displaystyle }& {\displaystyle }+M_{\mathrm{N}}\left[(L_{\mathrm{T}}^{\mathrm{2}}+z_{\text{NG}}{)}^{\mathrm{2}}+(q_{\mathrm{T}}+x_{\text{NG}}{)}^{\mathrm{2}}\right.\\ \text{(E30)}& {\displaystyle }& {\displaystyle }\left.+\mathrm{2}{\mathit{\nu }}_{\mathrm{y}}{q}_{\mathrm{T}}(z_{\text{NG}}{q}_{\mathrm{T}}-L_{\mathrm{T}}{x}_{\text{NG}})\right],{\displaystyle }{M}_{\mathrm{23}}=& {\displaystyle }{C}_{rT\mathrm{1}y}+\left[J_{\text{y,N}}+M_{\mathrm{N}}(x_{\text{NG}}^{\mathrm{2}}+z_{\text{NG}}^{\mathrm{2}}\right.\\ {\displaystyle }& {\displaystyle }\left.+L_{\mathrm{T}}{z}_{\text{NG}}+{\mathit{\nu }}_{\mathrm{y}}{q}_{\mathrm{T}}(z_{\text{NG}}{q}_{\mathrm{T}}-L_{\mathrm{T}}{x}_{\text{NG}})\right]{\mathit{\nu }}_{\mathrm{y}}\\ \text{(E31)}& {\displaystyle }& {\displaystyle }+M_{\mathrm{N}}\left[L_{\mathrm{T}}+z_{\text{NG}}\right],{\displaystyle }{M}_{\mathrm{33}}=& {\displaystyle }{M}_{\mathrm{e}}+M_{\mathrm{N}}+\left[J_{\text{y,N}}+M_{\mathrm{N}}\left(x_{\text{NG}}^{\mathrm{2}}-\mathrm{2}{x}_{\text{NG}}{q}_{\mathrm{T}}\right.\right.\\ \text{(E32)}& {\displaystyle }& {\displaystyle }\left.\left.+z_{\text{NG}}^{\mathrm{2}}\right)\right]{\mathit{\nu }}_{\mathrm{y}}^{\mathrm{2}}+\mathrm{2}{M}_{\mathrm{N}}{\mathit{\nu }}_{\mathrm{y}}{z}_{\text{NG}}.\end{array}$$

The elements of the forcing vector are

$$\begin{array}{ll}{\displaystyle }{f}_{\mathrm{1}}=& {\displaystyle }{f}_{\mathrm{H}}+\left[M_{\mathrm{F}}{z}_{\text{FG}}-M_{\mathrm{d}z}+M_{\mathrm{N}}(L_{\mathrm{T}}\right.\\ {\displaystyle }& {\displaystyle }\left.+z_{\text{NG}}-{\mathit{\nu }}_{\mathrm{y}}{x}_{\text{NG}}{q}_{\mathrm{T}})\right]{\mathit{\varphi }}_{\mathrm{y}}{\dot{\mathit{\varphi }}}_{\mathrm{y}}^{\mathrm{2}}\\ \text{(E33)}& {\displaystyle }& {\displaystyle }+M_{\mathrm{N}}\left[q_{\mathrm{T}}+x_{\text{NG}}+{\mathit{\nu }}_{\mathrm{y}}{z}_{\text{NG}}{q}_{\mathrm{T}}\right]{\dot{\mathit{\varphi }}}_{\mathrm{y}}^{\mathrm{2}}{\displaystyle }& {\displaystyle }+\left[\mathrm{2}{C}_{\text{tx}}+M_{\mathrm{N}}(\mathrm{1}+{\mathit{\nu }}_{\mathrm{y}}{z}_{\text{NG}}-{\mathit{\nu }}_{\mathrm{y}}^{\mathrm{2}}{x}_{\text{NG}}{q}_{\mathrm{T}})\right]\\ \text{(E34)}& {\displaystyle }& {\displaystyle }\times {\mathit{\varphi }}_{\mathrm{y}}{\dot{\mathit{\varphi }}}_{\mathrm{y}}{\dot{q}}_{\mathrm{T}}+M_{\mathrm{N}}{\mathit{\nu }}_{\mathrm{y}}\left[x_{\text{NG}}+{\mathit{\nu }}_{\mathrm{y}}{z}_{\text{NG}}{q}_{\mathrm{T}}\right]{\dot{\mathit{\varphi }}}_{\mathrm{y}}{\dot{q}}_{\mathrm{T}}\text{(E35)}& {\displaystyle }& {\displaystyle }+M_{\mathrm{N}}{\mathit{\nu }}_{\mathrm{y}}^{\mathrm{2}}\left[x_{\text{NG}}+z_{\text{NG}}{\mathit{\varphi }}_{\mathrm{y}}\right]{\dot{q}}_{\mathrm{T}}^{\mathrm{2}}\text{(E36)}& {\displaystyle }& {\displaystyle }+f_{\mathrm{a}}\left[\mathrm{1}-{\mathit{\theta }}_{\mathrm{t}}{\mathit{\nu }}_{\mathrm{y}}{q}_{\mathrm{T}}-{\mathit{\nu }}_{\mathrm{y}}{\mathit{\varphi }}_{\mathrm{y}}{q}_{\mathrm{T}}\right],\text{(E37)}& {\displaystyle }{f}_{\mathrm{2}}=& {\displaystyle }{\mathit{\tau }}_{\mathrm{H}}+M_{\mathrm{N}}\left[{\mathit{\nu }}_{\mathrm{y}}^{\mathrm{2}}(L_{\mathrm{T}}{x}_{\text{NG}}-z_{\text{NG}}{q}_{\mathrm{T}})\right]{\dot{q}}_{\mathrm{T}}^{\mathrm{2}}{\displaystyle }& {\displaystyle }-\mathrm{2}{M}_{\mathrm{N}}\left[q_{\mathrm{T}}+x_{\text{NG}}+{\mathit{\nu }}_{\mathrm{y}}(\mathrm{2}{z}_{\text{NG}}{q}_{\mathrm{T}}-L_{\mathrm{T}}{x}_{\text{NG}})\right.\\ \text{(E38)}& {\displaystyle }& {\displaystyle }\left.-{\mathit{\nu }}_{\mathrm{y}}^{\mathrm{2}}{q}_{\mathrm{T}}(L_{\mathrm{T}}{z}_{\text{NG}}+x_{\text{NG}}{q}_{\mathrm{T}})\right]{\dot{\mathit{\varphi }}}_{\mathrm{y}}{\dot{q}}_{\mathrm{T}}{\displaystyle }& {\displaystyle }+g\left[M_{\mathrm{F}}{z}_{\text{FG}}{\mathit{\varphi }}_{\mathrm{y}}-M_{\mathrm{d}z}{\mathit{\varphi }}_{\mathrm{y}}+M_{\mathrm{N}}\left\{(L_{\mathrm{T}}+z_{\text{NG}}\right.\right.\\ \text{(E39)}& {\displaystyle }& {\displaystyle }\left.\left.-{\mathit{\nu }}_{\mathrm{y}}{x}_{\text{NG}}{q}_{\mathrm{T}}){\mathit{\varphi }}_{\mathrm{y}}+q_{\mathrm{T}}+x_{\text{NG}}+{\mathit{\nu }}_{\mathrm{y}}{z}_{\text{NG}}{q}_{\mathrm{T}}\right\}\right]{\displaystyle }& {\displaystyle }+f_{\mathrm{a}}\left[L_{\mathrm{T}}+z_{\text{NR}}+{\mathit{\theta }}_{\mathrm{t}}{x}_{\text{NR}}+{\mathit{\theta }}_{\mathrm{t}}{q}_{\mathrm{T}}+{\mathit{\nu }}_{\mathrm{y}}{q}_{\mathrm{T}}^{\mathrm{2}}\right.\\ \text{(E40)}& {\displaystyle }& {\displaystyle }\left.-L_{\mathrm{T}}{\mathit{\theta }}_{\mathrm{t}}{\mathit{\nu }}_{\mathrm{y}}{q}_{\mathrm{T}}\right],\text{(E41)}& {\displaystyle }{f}_{\mathrm{3}}=& {\displaystyle }{f}_{\mathrm{e}}-D_{\mathrm{e}}{\dot{q}}_{\mathrm{T}}-K_{\mathrm{e}}{q}_{\mathrm{T}}{\displaystyle }& {\displaystyle }+M_{\mathrm{N}}\left[q_{\mathrm{T}}+x_{\text{NG}}+{\mathit{\nu }}_{\mathrm{y}}(\mathrm{2}{z}_{\text{NG}}{q}_{\mathrm{T}}-L_{\mathrm{T}}{x}_{\text{NG}})\right.\\ \text{(E42)}& {\displaystyle }& {\displaystyle }\left.-{\mathit{\nu }}_{\mathrm{y}}^{\mathrm{2}}{q}_{\mathrm{T}}(L_{\mathrm{T}}{z}_{\text{NG}}+x_{\text{NG}}{q}_{\mathrm{T}})\right]{\dot{\mathit{\varphi }}}_{\mathrm{y}}^{\mathrm{2}}\text{(E43)}& {\displaystyle }& {\displaystyle }+M_{\mathrm{N}}{\mathit{\nu }}_{\mathrm{y}}^{\mathrm{2}}{x}_{\text{NG}}{\dot{q}}_{\mathrm{T}}^{\mathrm{2}}{\displaystyle }& {\displaystyle }+g\left[C_{tT\mathrm{1}x}{\mathit{\varphi }}_{\mathrm{y}}+M_{\mathrm{N}}\left({\mathit{\nu }}_{\mathrm{y}}{x}_{\text{NG}}+{\mathit{\nu }}_{\mathrm{y}}^{\mathrm{2}}{z}_{\text{NG}}{q}_{\mathrm{T}}\right.\right.\\ \text{(E44)}& {\displaystyle }& {\displaystyle }\left.\left.-{\mathit{\nu }}_{\mathrm{y}}^{\mathrm{2}}{x}_{\text{NG}}{\mathit{\varphi }}_{\mathrm{y}}{q}_{\mathrm{T}}+{\mathit{\nu }}_{\mathrm{y}}{z}_{\text{NG}}{\mathit{\varphi }}_{\mathrm{y}}+{\mathit{\varphi }}_{\mathrm{y}}\right)\right]\text{(E45)}& {\displaystyle }& {\displaystyle }+f_{\mathrm{a}}\left[\mathrm{1}+{\mathit{\theta }}_{\mathrm{t}}{\mathit{\nu }}_{\mathrm{y}}{x}_{\text{NR}}-{\mathit{\theta }}_{\mathrm{t}}{\mathit{\nu }}_{\mathrm{y}}{q}_{\mathrm{T}}+{\mathit{\nu }}_{\mathrm{y}}{z}_{\text{NR}}\right].\end{array}$$