A symbolic framework for flexible multibody systems applied to horizontal axis wind turbines

Abstract. The article presents a symbolic framework that is used to obtain the linear and non-linear equations of motion 1 of a multibody system including rigid and flexible bodies. Our approach is based on Kane’s method and a nonlinear shape 2 function representation for flexible bodies. The method yields compact symbolic equations of motion with implicit account of 3 the constraints. The general and automatic framework facilitate the creation and manipulation of models with various levels 4 of fidelity. The symbolic treatment allows for the obtention of analytical gradients and linearized equations of motion. The 5 linear and non-linear equations can be exported to Python code or dedicated software. The application are multiple such as: 6 time-domain simulation, stability analyses, frequency domain analyses, advanced controller design, state observers, digital 7 twins, etc. In this paper, we describe the method we used to systematically generate the equations of motion of multibody 8 systems. We apply the framework to generate illustrative onshore and offshore wind turbine models. We compare our results 9 with OpenFAST simulations and discuss the advantages and limitations of the method. A Python implementation is provided 10 as an opensource project. 11

acceleration of the point are given by (see e.g. Shabana (2013)): 48 r P = r i + s P = r i + s P0 + u P (1) where: r i , v i and a i are the position, velocity and acceleration of the origin of the body; s P0 is the initial (undeformed) 52 position vector of point P with respect to the body origin; u P is the elastic displacement of the point (0 for rigid bodies); 53 ω i is the rotational velocity of the body with respect to the inertial frame; (˙) and (˙) i refer to time derivatives in the inertial 54 f r + f * r = 0, r = 1 . . . n q (4) 62 where f * r is associated with inertial loads and f r is associated with external loads, and these components are obtained for each 63 generalized coordinates. The components are obtained as a superposition of contributions from each body: The terms f ri and f * ri can be obtained for each body individually and assembled at the end to form the final system of equa-66 tions. We will present in subsection 2.3 and subsection 2.4 how these terms are defined for rigid bodies and flexible bodies 67 respectively. 68 2.3 Rigid bodies 69 We assume that body i is a rigid body and proceed to define the terms f ri and f * ri . The inertial force, f * i , and inertial torque, 70 τ * i , acting on the body are: where m i is the mass of the body, a G,i is the acceleration of its center of mass with respect to the inertial frame, I G,i is the 73 inertial tensor of the body expressed at its center of mass. Equation 6 is a vectorial relationship, it may therefore be evaluated 74 in any coordinate system. The component f * ri is defined as: 76 with 77 J v,ri = ∂v G,i ∂q r , J ω,ri = ∂ω i ∂q r (8) 78 where v G,i is the velocity of the body masscenter with respect to the inertial frame. The partial velocities, or Jacobians, J v 79 and J ω , are key variables of the Kane's method. They project the physical coordinates into the generalized coordinates (q), 80 inherently accounting for the kinematic constraints between bodies. In numerical implementations, the Jacobians are typically 81 J v,ri = ∂v O,i ∂q r , J ω,ri = ∂ω i ∂q r , J e,ri = ∂q e,i ∂q r (20) 120 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 In various applications, a linear time invariant approximation of the system is desired. Such approximation is obtained at an operating point, noted with the subscript 0, which is a solution of the non-linear equations of motion, viz: 136 e(q 0 ,q 0 ,q 0 , u 0 , t) = 0 (23) 137 The linearized equations about this operating point are obtained using a Taylor-Series expansion: 138 M 0 (q 0 )δq + C 0 (q 0 ,q 0 , u 0 )δq + K 0 (q 0 ,q 0 , u 0 )δq = Q 0 (q 0 ,q 0 , u 0 )δu  3 Implementation into a symbolic framework 145 In this section we discuss a Python open-source symbolic calculation framework that implements the equations given in sec-146 tion 2. A Maxima implementation from the same authors is also available Geisler (2021).

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

156
For the FlexibleBody class, we followed the formalism of Wallrapp (1994), and implemented Taylor expansions for all body is assumed to occur at one extremity of the flexible body. Some knowledge of SymPy mechanics are still required to use this layer.
3) The third level consists of template models such as generic onshore or offshore wind turbine models. Degrees of 166 freedom are easily turned on and off for these models depending on the level of fidelity asked by the user and generic external 167 forces can be implemented or declared as external inputs. The non-linear and linear equations of motion can be exported to 168 latex and python-ready scripts for various applications (see subsection 5.1). As little as three lines of code are required by 169 the user to perform the full step from derivation of the equations, optional linearization, and exportation. To obtain numerical 170 results from the exported python code, the user needs to provide the arrays with the degrees of freedom values q andq, their 171 initial conditions,a dictionary with inputs (u) that are function of time, and a dictionary of parameters (p) containing all the 172 numerical constants such as mass, acceleration of gravity, etc. We provide tools to readily time-integrate the generated model 173 using numerical values (including initial values) from a set of OpenFAST input files.

174
The source code of YAMS is available on github as a subpackage of the Wind Energy LIBrary, WELIB (Branlard, 2021   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 197 to familiarize the reader with the key concepts of the shape function approach given in Appendix B. A sketch of the system 198 is given in Figure 1. We start by modelling the blade using a single shape function, assumed to be directed along the x axis 199 ("flapwise"): Φ 1 = Φx, where the hat notation indicates the unit vector in the x direction. The undeflected blade is directed The generalized force is obtained from Equation B3: The geometrical stiffness contribution of the axial load is obtained from Equation B5 as: The axial stiffness, K g , is positive and increases with the square of the rotational speed. This restoring effect is referred to as 215 "centrifugal stiffening". The natural frequency of the blade will increase with the rotational speed as follows: where k Ω is referred to as the rise factor, or Southwell coefficient, and in our approximation, it is found to be constant: The coefficient provides the variation of the blade frequency with rotational speed, which is something 219 that is observed on a Campbell diagram when performing stability analyses. In general, the mode shapes of the blade will also 220 change as function of the rotational speed, and different shape functions should preferably be used for simulations at different is negligible as long as the rotational speed is small compared to the natural frequency (e.g., (Ω/ω) 2 5, see Bielawa (2006)), 224 which is the case for wind energy applications.

225
The treatment for a shape function in the edgewise direction is similar, using Φ 2 = Φ 2θ . In this case, the centrifugal force 226 also has a component in the tangential direction equal to p θ,centri (r) = −Ω 2 u θ (r)dm(r), with u θ = Φ 2 q. This leads to a gen-227 eralized force equal to L 0 p θ,centri drΦ 2 = −Ω 2 M e q, or equivalently, to a stiffness term: K ω = −Ω 2 M e . It can be verified that 228 this generalized force corresponds to the contribution O e,11 ω 2 x , from k ω,e , given in Equation A10. For an edgewise mode, the 229 frequency therefore evolves as: with k Ω = K g (Ω)/M e /Ω 2 and with K g computed using Equation 30. 232 We apply the method to the NREL-5MW turbine using the blade properties and shape functions provided in the ElastoDyn  We consider a system of three bodies: tower (or support structure), nacelle and rotor. The system represents an onshore wind 247 turbine or a fixed-bottom offshore wind turbine. A sketch of the system is given in Figure 3. The nacelle and rotor blades are 248 rigid bodies, whereas the tower is flexible and represented by one shape function 2 in the fore-aft direction, noted Φ 1 = Φ 1x .

249
For hand calculations and as a first approximation, the first mode shape of a massless beam with a top-mass may be used: 250 Φ 1 (z) = 1 − cos(zπ/L/2). Increased accuracy is obtained when the shape function matches the actual first tower fore-aft e ω e where m(z) and 261 EI(z) are the mass per length and bending stiffness of the tower, and ω e and ζ are the frequency and damping ratio associated 262 with the shape function (assuming the shape function is close to a mode shape). The geometrical softening of the tower due to 263 the tower top mass (K gt ) and its self-weight (K gs ) is obtained using Equation B5, as K g = K gt + K gs , with : The shape function frequency is obtained as: The application of the symbolic framework leads to the following equations of motion (rearranged for interpretability): 274 obtained by taking the gradient of the forcing with respect to q, and using a small angle approximation for ν y to the second 294 order: L T (2a 2 z + 3a 3 z 2 + 4a 4 z 3 + 5a 5 z 4 + 6a 6 z 5 )/(a 2 + a 3 + a 4 + a 5 + a 6 )   We consider the same system as the one presented in subsection 4.4 but the tower is now represented by one shape function in 312 both the fore-aft and side-side directions, Φ 1 = Φ 1x and Φ 2 = Φ 2ŷ . The degrees of freedom are q = (q 1 , q 2 , ψ), where q 1 and 313 q 2 are the generalized (elastic) coordinates in the fore-aft and side-side directions respectively, and ψ is the rotor azimuth. A 314 sketch of the system is given in Figure 3.  in azimuthal speed, resulting from the coupling between the gyroscopic loads and the tower bending, is well captured.

333
In this example, we demonstrate the applicability of the method for a floating wind turbine. We model the turbine using 3 335 bodies: rigid floater, flexible tower, rigid rotor-nacelle-assembly (labelled "N"). The degrees of freedom selected are: q = 336 (x, φ, q T ), where x is the floater surge, φ, the floater pitch, and q T , the coordinate associated with a selected fore-aft shape 337 function. A sketch of the model is given in Figure 6. The notations are similar to the ones presented in subsection 4.5. Lumped of θ t and φ y , and a second order approximation for ν y .

342
We performed a numerical simulation of the model generated by YAMS and compared it with OpenFAST, for a case with 343 gravitational loads only, starting with x = 0 m, φ = 2 deg and q T = 1 m. The results are presented in Figure 7. We observe 344 again that the results from the two models correlate to a high degree. 345 We also compared the linearized version of both models. The symbolic framework can generate the linearized mass, stiffness

349
The 4% error in the pitch frequency appears reasonable in view of the approximations used.

383
-The expression of the displacement field u in terms of a superposition of shape function is typically done using a first 384 or second order expansion. We discuss these formulations in Appendix C. Our framework supports both approaches: it 385 is the responsibility of the user to provide the terms and values of the expansion when numerical evaluation is to occur.

386
The advantage of the second-order expansion is that some geometrical non-linearities are directly accounted for by the 387 method, whereas these non-linearities need to be introduced "manually" in the first-order expansion approach, as done  In spite of the advantages listed in subsection 5.2, the symbolic procedure presented in this work has some potential limitations. 404 We are identifying two in this section. First, constraints and closed loops have currently not been added to the framework. The The different terms of the mass matrix are obtained as follows: The integrals are understood as volume integrals over the volume of the body (for beams they reduce to line integrals). The coefficients (Wallrapp (1994)).
The body elastic forces are given by: where K e and D e are the elastic stiffness and damping matrix, and k σ represent stress stiffening terms. The elastic damping 464 forces are often given as stiffness proportional damping. For more details, see Wallrapp (1994), and for more example of with 465 elastic beams see Branlard (2019). The external loads can be assumed to consists of distributed volume forces p (in practice 466 they are mostly surface forces of line forces) and a gravitational acceleration field g . The components of the external loads in 467 Equation A1 are then obtained by integration over the whole body: Appendix B: Application of the shape function approach to an isolated beam    (2019)), the Jacobians are sometimes stacked into a matricial form: Some implementation choices are needed depending if these matrices are expressed in the global frame or a body frame. The

513
Jacobian matrices are referred to as velocity transformation matrix, and the link between formulations in global and local 514 coordinates is given in Branlard (2019). In the same reference, recursive relationships are given for tree-like assembly of 515 bodies, to help express the Jacobian matrices of each body recursively, based on the matrices of the parent body. It is also noted 516 that the quadratic velocity terms, k ω , can be obtained using the time derivative of the Jacobian matrix.

517
In this article, we have not explicitly written the rotational impact of the shape functions. In most applications, bodies are 518 connected at their extremities and the deflection slope at a body extremity will induce a rotation of the subsequent body. The 519 deflection slope can be obtained form the knowledge of the shape functions. This is readily accounted for by introducing time-520 varying rotation matrix between bodies, and this is the approach used in our symbolic framework. A formalism of rotations 521 of bodies connected at their extremities is given in Branlard (2019). A more general formulation, introducing shape function 522 rotations Ψ, is given in (Wallrapp, 1994;Schwertassek and Wallrapp, 1999;Lemmer, 2018). In such formulation, the linear 523 rotation field is obtained as I + Ψq, where I is the identity matrix.

524
The order of expansion of the displacement field leads to alternative formulations. In Shabana (2013) and Branlard (2019) a 525 first order expansion is used: u = j [Φ 0 j ]q e,j In the work of Wallrapp a second order expansion is used: u = j Φ 0 j + 1 2 k Φ 1 jk q e,k q e,j

526
In both formulations, the equations of motion given in Appendix A lead to shape-integral expansions of the following form: where T is a dummy variable standing for M θ,θ , C t , C r , G r , G e , or O e . The "0" and "1" terms are stored using a "Taylor" 529 object-oriented class in the SID format defined by Wallrapp. The subtlety lays in the fact that the "1" terms will be different 530 if the displacement is developed using a first order expansion or a second order expansion. Some terms involving Φ 1 jk will be present in the latter case. The reader is referred to Wallrapp (1993) for a full description of the Taylor-expanded terms. Setting Φ 1 jk = 0 in these expressions will lead to the 1st-order shape integral approach of Shabana. The nacelle kinematics (at its center of mass) are: a G,N =qt x + (−ν 2 y x N Gq 2 + ν y z N Gq )n x + (−ν 2 y z N Gq 2 − ν y x N Gq )n z (D4)

543
The Jacobians with respect to q are: