Implementation of an advanced beam model in BHawC

The implementation of a new structural beam model in the in-house aeroelastic tool BHawC is presented. The model is based on a recently published equilibrium based formulation of a beam stiffness matrix as well as a novel stiffness proportional damping formulation. The performance of the new beam implementation is assessed via modelling a commercial wind turbine. Four analyses are presented to illustrates the increased accuracy in recovered section forces and representation of the dynamic behaviour gained from the new structural model.


Introduction
Design bases for wind turbines require a significant number of load simulations, for offshore projects this ranges in the tens of thousands. Reducing the time each simulation can both reduce the overall computational load and reduce the time the engineer has to wait in order to ensure that the results were as expected and possibly perform another iteration. One way to reduce the time a simulation takes is to reduce the order of the structural model. This paper presents an implementation of equilibrium beam elements which allows using fewer elements while retaining the accuracy previously obtained by using a higher number of simpler prismatic beam elements. Key factors in allowing the discretization of a substructure with fewer elements are that section force results can be extracted anywhere inside an element with great accuracy and that the stiffness and damping matrices consistently account for lengthwise variations in properties. The separation of the discretization of the model from the definition of the model has the additional practical benefit that additional results can be extracted at a later point without introducing small variations in the results, which would otherwise happen if nodes were to be added or moved.
The new beam element formulation presented in this paper is composed of two parts, namely the stiffness matrix and section force recovery formulation presented by Krenk and Couturier [2], and a novel stiffness proportional damping formulation. The underlying theory is based on a complementary-energy formulation for a two-node straight nonhomogeneous anisotropic elastic beam element which makes use of equilibrium internal-force distributions without any need for shape functions. This theory has the advantage of being able to accept directly 6 by 6 cross-section stiffness matrices and account for arbitrary lengthwise variation in properties. This makes it well suited for the design of blades with deformation mode coupling such as bend-twist coupling. Furthermore, the section force recovery approach allows obtaining the true value of any of the six section forces anywhere inside an element.
It is convenient when modelling the damping of non-homogeneous structures, such as fiber reinforced wind turbine blades, to be able to control the level of damping of the different deformations independently. The well known Rayleigh damping provides a very convenient way of introducing damping as it only requires the tuning of two parameters which scale the entire mass and stiffness matrices, respectively. In a aeroelastic simulation codes where the degrees of freedom are absolute and defined in an inertial frame, mass-proportional damping introduces damping of rigid body motion of the elements which is not the intention. The use of only a single parameter on the stiffness matrix however prevents the tuning of different deformations. Control over the damping of each mode of the structure can be achieved using a modal damping approach. This approach however tends to increase calculation effort when running a nonlinear analysis as it reduces the sparsity of the damping matrix and requires performing several eigenvalue analyses.
The novel stiffness proportional damping formulation presented in this paper represents a natural extension of the equilibrium based stiffness matrix formulation. Similar to the work of Hansen [3], individual damping coefficients are used to control the flap, edge, torsion and extension damping. The damping coefficients in the present approach are applied to the crosssection flexibility matrix. Integration of the static field description over the element length provides the element damping matrix directly. The resulting damping matrix is invariant to the beam element coordinate system used and accounts for arbitrary lengthwise stiffness variations.
The equilibrium beam element formulation has been implemented in the in-house aeroelastic simulation code BHawC which has been used in Siemens Gamesa Renewable Energy for almost 15 years for producing loads for the design of wind turbines. The code is being continuously validated against measurements from prototype turbines. The structural model in BHawC is built around the co-rotational formulation providing geometric nonlinearity [1]. The substructures for the tower, shaft and blades are modelled using beam elements which are defined by continuous cross-section properties tabulated either in a reduced or in a 6 by 6 cross-section format.
The next section describes the theory behind the equilibrium beam element and defines the stiffness and damping matrices as well as the internal-force recovery method. Then the application of the theory is illustrated through four numerical examples using a model of a Siemens Gamesa wind turbine.

Equilibrium beam element
This section presents the complementary formulation for a linear two-node straight beam element which has recently been implemented in BHawC. The presentation of the theory is split in three sections. The first two sections summarize the derivation for the stiffness matrix, the equivalent nodal forces, and section force distribution recovery developed by Krenk and Couturier [2]. The third section presents the formulation for a new element stiffness proportional damping matrix which represents a natural extension of the complementary based stiffness matrix.  The theory is based on the analysis of a straight beam with length l, longitudinal coordinate x 3 , and cross-section coordinates x 1 and x 2 , as shown in figure 1. The beam static field can be characterized in terms of the distribution of the six generalized section forces and moments at each cross-section plane along x 3 which are grouped together in the force vector q( are two shear forces, and Q 3 (x 3 ) is the axial force. The components M 1 (x 3 ) and M 2 (x 3 ) are bending moments, and M 3 (x 3 ) is the torsion moment component. Unlike the kinematic field description of a beam using shape functions which introduces complications when treating non prismatic beams, the statics of a beam can be described by six simple equilibrium modes, namely the homogeneous state of extension, shear, bending, and torsion. Moreover, the distribution of the generalized internal forces and moments q h (x 3 ) for a beam with no external loads can be represented by the following relation where q 0 is the force vector representing the six generalized section forces at the middle, l/2, of the beam. The 6 by 6 interpolation matrix T(x 3 ) provides the spatial linear variation of the bending moments along x 3 .

Stiffness matrix
Linear stiffness matrix formulations using classic kinematics or equilibrium based approaches can be developed by integrating the cross-section strain energy over the element's length. In classic kinematic formulations, the strains, the curvatures and the cross-section stiffness properties along x 3 are required. Advanced cross-section analysis solvers [4,5,6] most commonly provide this information in terms of cross-section 6 by 6 stiffness matrices D. In the present equilibrium based formulation, the strain energy is calculated using a complementary energy approach based on integrating the static field q h (x 3 ) and the associated cross-section flexibility matrix C = D −1 , The strain energy on the right is expressed in terms of the mid-beam force vector q 0 via the static field described in (1) and the beam element flexibility matrix H associated with the six equilibrium modes defined by the integral The strain energy of a beam element with six equilibrium modes can also be written in terms of the element 12 by 12 stiffness matrix K and the 12 nodal displacements and rotations ordered in the vector u as From energy equivalence of (2) and (4) it can be shown that the beam element stiffness matrix K can be expressed in terms of the beam element flexibility matrix H and end point values of The resulting stiffness matrix formula directly accounts for varying cross-section properties, deformation mode coupling such as bend-twist coupling, and transverse shear deformation without the need to calculate shape functions.

Internal section force distribution and equivalent nodal load
The distribution of internal section forces q(x 3 ) in a beam with external loads is composed of two parts, namely the distribution of section forces corresponding to the homogeneous solution q h (x 3 ) presented in (1), and the distribution of forces which are in equilibrium with the applied external loadq( The internal force distributionq(x 3 ) is not unique and can be chosen to represent the equilibrium state of the beam element subject to the external load for any combination of support conditions at the nodes. The choice of a different boundary condition to obtainq(x 3 ) will change the homogeneous solution q h (x 3 ) without effecting the total internal section force distribution q(x 3 ).
The total potential energy of the beam element Π p expressed in terms of the internal force distribution q(x 3 ) can be expressed in a compact two part form with the strain energy of the beam element expressed in terms of the beam element stiffness matrix K, and the potential energy of the applied load in terms of the equivalent nodal forces r, The equivalent nodal load vector r is given by where end point values corresponding to the internal force distributionq( It is seen that the present energy consistent equivalent nodal load vector accounts for both arbitrary applied distributed load and arbitrary lengthwise variation of the cross-section properties. From the principle of stationary potential energy of the beam element, the expression for the homogeneous part of the internal force distribution q h (x 3 ) comes out as The total distribution of internal section forces q(x 3 ) follows from (6) by summing the homogeneous part of the force distribution from (10) and the internal force distributionq(x 3 ) calculated at the beginning of the simulation. Section forces calculated from the conventional finite element method are limited to constant and linear variations of section forces within an element as they only capture the homogeneous part of the total distribution of internal section forces.

Stiffness proportional damping matrix
It is convenient in practical applications of beam models to be able to control the level of damping of the different deformation modes independently. Hansen [3] developed a convenient approach to control the damping of torsion and bending of prismatic beams using a stiffness proportional damping matrix obtained by scaling the three moments of inertia independently. Extending this concept to non-uniform beams can be approached in two ways, namely applying scaling factors on the element stiffness matrix or at the cross-section level. For both approaches, the scaling needs to be applied in the principal coordinate system. This poses a certain difficulty when applying the scaling factors at the element level as the calculation of the element's principal axis is non-trivial for beams with varying cross-section properties such as structural pitch. In the present formulation, scaling factors are applied on the cross-section flexibility matrix C directly in order to obtain a damped 6 by 6 cross-section flexibility matrix C d , Contrary to the element stiffness matrix, expressing the cross-section flexibility matrix C in its principal coordinate can easily be done using standard mechanics of structures theory. The main diagonal of the scaling matrix S contains four scaling coefficients It can be seen that the damping coefficients b 1 and b 2 are used to scale the bending deformations, the coefficient b 3 scales the torsion deformation, and b 4 scales the extension deformation. In cases where the extension mode is of little importance, as is the case for current wind turbine blades, the extension scaling factor can be set to a constant value. This limits the number of damping variables that need to be found experimentally. At this point, the element stiffness proportional damping matrix calculation follows the stiffness matrix formulation summarized in section 2.1. Substituting the cross-section flexibility matrix C for the damped cross-section flexibility matrix C d in (3) gives the damped beam flexibility matrix The beam stiffness proportional damping matrix B follows from substituting the damped beam flexibility matrix H d in for the flexibility matrix H in (5), The resulting damping matrix is invariant to the beam element coordinate system and accounts for arbitrary lengthwise stiffness property variations including structural pitch. The matrix is guaranteed to remain semi-positive definite as the damping related scaling is introduced by pre and post multiplying the semi-positive definite flexibility matrix C by a diagonal matrix in (11) which does not change the eigenvalues of the matrix.

Analysis
This section illustrates the application of the present new structural model in BHawC to analyse commercial wind turbines. The first example compares recovered axial force in a tower using conventional finite element approach, the new equilibrium approach, and the analytical solution.
The second example shows the convergence of the modal frequencies and damping for the standard and equilibrium models. The third example shows the distribution of all six recovered section forces along the blade compared to the forces obtained using the conventional finite element nodal method. The last example compares the dynamic behaviour of a full wind turbine using the different methods.

Tower section force recovery
The first example concerns the analysis of a 118 m wind turbine steel tower. The tower is discretized using 13 beam elements. The distribution of the axial force along the tower at standstill is shown in figure 2. The distributions are normalized with respect to the axial force at the tower bottom. It can be seen that the recovered section force distribution using the equilibrium method agrees perfectly with the analytical solution obtained by integrating the gravity load from the mass distribution directly. Conversely, the piecewise constant force distribution obtained from the conventional finite element method at the nodes fails to accurately capture the true continuous axial force distribution. The conventional finite element force recovery approach under predicts the tower bottom compressive force by 3 %. A workaround for this limitation is to include a short element around the point where one would like to extract the force precisely, as illustrated at the tower top. This is not necessary with the section force recovery of the equilibrium formulation which prevents increasing the number of degrees of freedom to solve.

Convergence of blade modal frequencies and damping
The next example illustrates the convergence properties of the standard and the equilibrium beam models for a 75-m-long wind turbine blade currently manufactured by Siemens Gamesa Renewable Energy. Figure 3 shows the difference in modal frequency and damping for different numbers of equal-length beam elements compared to the reference simulations with 128 elements for each model. The three modes shown, namely the first flapwise bending, the first edgewise bending, and the first torsional modes, are tuned to the same logarithmic decrement with both models using 128 elements. The scaling coefficients b 1 , b 2 , b 3 and b 4 of (12) thus obtained are kept the same for the other discretizations. The frequencies for these three modes obtained using 128 elements agree within 0.1 % between the two models, which is considered to be a converged level. The error in the frequency for the first flapwise and edgewise modes using the equilibrium model is within 0.1 % even when using four elements whereas it is approximately within this converged level only when using 32 or more elements with the standard model. The frequency of the first torsional mode has a higher sensitivity to the element discretization because it also contains some amount of higher-order flapwise and edgewise bending shapes which cannot be captured with too few elements. Using the equilibrium model, the damping of the three modes are within 1 % error when using 16 elements whereas 64 elements are required with the standard 7 1234567890 ''""  It should be noted that for a complex structure such as the current blade more accurate values of modal frequencies and damping for both the standard and the equilibrium models can be obtained by optimizing the beam node locations according to the most relevant mode shapes. However, the tendency of faster convergence for the equilibrium model is still present when using optimized node locations.
The convergence rate of the dynamic modes of the equilibrium model allowed to reduce by half the number of beam elements used to model a full turbine without affecting the loads. Running with such a smaller model allowed to reduce approximately also by half the computational time of a simulation. Note that the computational gain depends on the profile of the program, i.e. which tasks are taking the most time. In the current implementation, all the matrices associated with the new beam formulation are solved only once during the initialization phase thereby having no noticeable effect on the computational time.

Blade section force recovery
This example concerns the same blade as in the previous example. The blade is discretized using eight beam elements. The nodes are positioned along the elastic axis and their lengthwise distribution is optimized to minimize the error of the first four blade natural frequencies and the first torsion frequency.
The internal distribution of all six section forces recovered using the present equilibrium method and the standard finite element approach are shown in figure 4. The section forces are shown in the blade root coordinate system and are normalized with respect to the largest section forces calculated using the equilibrium method. This contemporaneous view of the section forces on a blade is an example of a loading situation from a dynamic simulation. As with the tower example, the piecewise constant nature of the shear and axial force distributions obtained from conventional finite element method under or over predicts the internal shear and axial forces. Even the piecewise linear variation of the bending moments under predicts the flap and edge root peak moments by 2 % compared to the moments recovered using the present equilibrium approach. Using a finer discretization would enable conventional finite element based section forces to converge towards equilibrium based recovered forces however at the cost of increasing computational time.

Dynamic Analysis
This final example of the present paper concerns the dynamic analysis of a commercial Siemens Gamesa wind turbine. The blade and tower used in the present analysis is the 75-m-long blade and 118-m high tower studied in the previous examples. The turbine is subjected to gust event with normal shut-down. In the simulation, the gust event starts at t = 20 sec and the turbine shuts down at t = 26 sec. Tree different models are used to represent the turbine. In the first model, 16 standard beam elements are used to discretize each blade. The second model is identical to the previous one but with half the number of elements per blade. The third model uses eight equilibrium beam elements to model each blade. Time series of the blade root flap moments and tower bottom fore-aft moments obtained using all three models are shown in figure 5. The results are normalized with respect to the largest moments calculated using the 16 element model. It can be seen that the moments obtained using the model with equilibrium based beam elements agrees well with the moments from the 16 element model. Conversely, the moments from the model with a reduced standard beam element mesh has a considerably different dynamic response to the gust and also yields a 5 % smaller peak flap moment. Using only eight standard beam elements changes the turbine 9 1234567890 ''""  properties leading to a different controller response to the gust which resulted in different loads. This illustrates a potential pitfall of reducing the mesh size when using standard beam elements.
Because of the higher model convergence rate, this issue would only occur with equilibrium beam elements when trying to use an even smaller mesh size.

Conclusions
A theory for a non-homogeneous two-node straight beam element previously published has been extended with the development of a new element stiffness proportional damping matrix. In the new formulation, individual damping coefficients applied to the cross-section flexibility matrices allow to control the flap, edge, torsion, and extension damping independently. The approach accepts six by six stiffness matrices and naturally accounts for any length-wise variations in stiffness properties including tapering, twist, and deformation mode coupling such as bendtwist coupling. The beam model has been implemented in the in-house aeroelastic code used at Siemens Gamesa called BHawC. The method shows increased accuracy modelling static and dynamic behaviours, illustrated by the analysis of a commercial wind turbine. The results highlight that the use of the new equilibrium based beam element can decrease the number of elements required to achieve model convergence thereby allowing substantial computational cost reduction.