Integration of electromagnetic finite element models in a multibody simulation to evaluate vibrations in direct-drive generators

This paper introduces a novel electromechanical model for calculating electromagnetic excited structural vibrations and structure borne acoustics for gearless wind turbines. Therefore, the wind turbine model structure is explained and a drivetrain model is derived to investigate the drivetrain decoupled from the aerodynamic excitations. The drivetrain model is fed with results from an electromagnetic finite element model of the generator considering air gap width changes and the wind turbine torque and speed characteristics. Furthermore, an exemplary ramp-up of the drivetrain is simulated. It can be seen, that generator structure oscillations are excited during certain rotational speeds, which may be relevant for the acoustic behavior of the turbine.

air gap forces excite the generators structural components, which causes vibration. These vibrations are even enhanced due to the change of airgap widths due to wind loads. Vibrations may cause excessive wear in drivetrain components, which may lead to plant downtime due to component failures [3][4][5]. Drivetrain vibrations also cause acoustic noise emission [6,7], which restricts the erection of new plants [8,9]. Therefore, a coupled model of the electromagnetic excitation and structural components of the turbine is needed to predict vibrations caused by the generator especially for direct-drive generators.
This paper introduces a multiphysical simulation approach for direct-drive generators, which integrates forces from the electromagnetic finite element solution into

Wind turbine simulation model
To model structural vibrations of a gearless wind turbine, considering wind loads and electromagnetic loads from the generator and the turbine power curve, a transient multibody model approach is chosen. The model results are the deformation and acceleration of different eigenmodes, which can be extracted and analyzed (see Sect. 5). These results can then be used as input data for an air borne acoustic simulation to calculate the sound emission of the turbine according to the field measurement described in the IEC61400-11, or to further optimize the power curve behavior or structural and electromagnetic components of the turbine [10].
The model approach for the wind turbine utilizes the multibody simulation environment Simpack (see Fig. 1). In this model, the turbines structural components are mod- eled as elastic bodies by means of modal reduction with ABAQUS [11]. The major sound emitting components are the rotor blades, the tower, the hub and nacelle cover. These components similarly are formulated as flexible bodies. As these components contain a large amount of eigenmodes the amount of modes are compressed to enhance the computational efficiency of the overall turbine model. For this purpose, a new method of modal compression was developed to compress the mode set of the flexible bodies as part of the ongoing research project. With the utilization of flexible bodies along with models for the main, pitch and yaw bearing in sufficient modeling depth the transfer of structureborne sound can be evaluated from the excitation locations to the sound emitting surfaces. Since the mode shapes of the sound emitting surfaces are well known from the modal reduction, these can be combined with the modal participation factors of the structures to regain the surface accelerations without direct measurement of all surface nodes which reduces the memory amount substantially. The sur-  face accelerations are then used to determine the emitted airborne sound [12]. The aerodynamic excitations are included in the model by the AeroDyn force element, which can impose the wind loads from certain wind profiles to the rotor blades. Main bearings, yaw bearing and pitch bearings are also included as force elements into the system model [13].
The generators air gap flux density is precalculated with the finite element tool iMOOSE (see Sect. 3). By using this air gap flux density the excitation forces of the generator are determined, including harmonics, saturation effects and locally varying air gap widths. These excitation forces are included by look-up tables. To impose the electromagnetic forces to the generator structure, a novel force element (EM-Force in Fig. 1) has been developed and is introduced in Sect. 4. This force element splits the air gap surfaces into subsurfaces and imposes the local forces depending on the rotors angular position and the local air gap widths.

Theory
In this chapter, the finite element-based calculation of electromagnetic force excitations is explained and derived.
As illustrated in Fig. 1, the electromagnetic modeling is embedded in the whole simulation tool chain in the form of look-up tables (LUT). These tables collect a priori-calculated electromagnetic forces in cylindrical coordinates (radial, tangential and axial) parametrized by 2D-finite element simulations using iMOOSE [14]. Assuming that the local field solution of one pole pitch geometry is not influenced by the adjacent poles, which holds true for the regarded multipole topologies, a respective partial model can be used. The basis for the calculation is the A-formulation of the magneto static problem by the magnetic vector potential A that is solved numerically for each node of the mesh with and where r is the nabla-operator, E B is the magnetic flux density, E J source is the electrical current density of the generator coils and is the material dependent magnetic reluctivity.
The last-named quantity shows for the magnetic flux leading components, that are the electrical steel stacks of the rotor and stator respectively, a non-linear behavior due to magnetic saturation. This leads in turn to a change in the spatial and temporal distribution of the magnetic air gap flux and therefore the exciting force densities at the material transition. The force density σ at the time instant t can be described in the Fourier-decomposed formulation where k is the spatial ordinal wave number and ! l is the angular frequency that is proportional to the rotational speed of the turbine in case of the regarded direct-driven synchronous generators. The orders k and l are the spatial and temporal orders respectively, k is the phase angle of the k-th spatial order.
Based on the formulation in (3) a mathematical description for the sum forces can be derived: F sum; i is the force exerted by the i-th stator tooth or rotor pole slice of length l slice and width b of N in total respectively. d is the bore diameter. The orientation of the forces can be radial and tangential, which physical origin can be derived from the simplified Maxwell stress tensor, correlating it with the magnetic air gap flux density, in case of 2D-models.

Modeling and results
In [15] it is shown that the electromagnetic forces depend beside the operating point on the air gap width. The sensitivity towards magnetic saturation along the characteristic operating curve of the wind turbine is shown at no-load conditions and the LUT-based approach is motivated to consider the non-linear effects described in the preceding section. Moreover, the axial varying electromagnetic field and force conditions are presented due to axial deformations of the air gap, as e.g. tilting of the generator rotor. The composition of the local forces by spatial distributed markers of the force element with air gap-recalculation is shown and are broadly presented in the following chapter.
In this work, the analysis of axial effects is extended by the impact of the rotor pole shoe skewing. Moreover, the magnetic forces are evaluated for electric loaded-conditions as the basis for the analysis by the force element in additionally deformed condition.
In summary, the following dimensions are consequently considered for the LUT: The respective sum forces are derived for both electromagnetic active parts of the generator, stator and rotor, based on the nodal forces of the finite element meshes using the eggshell method [16].
At first, the results for the modeling of the skewed pole shoes based on the multi-slice approach [17] are presented. In Fig. 2, the radial and tangential force exerted by a rotor pole divided in 5 slices depending on the rotation angle related to the pole pitch p is shown.
The values are normalized to the mean value of the middle slice no. 3. Both force components are impacted by skewing as the waveforms and particularly the mean value change. In case of the tangential force, the torque varies consequently depending on the axial direction. In conclusion, the skewing effect can be considered efficiently by embedding the results for the different force element slice markers in the LUTs. In the second and final part of this section, results for each one operating point in partial (P WT 0.1P N ) and in full load of the wind turbine are presented. The normalized radial and tangential components of the stator tooth forces (exerted by the middle slice) depending on the angular rotor position for each load condition are illustrated in Fig. 3.
Both curves are normalized to the absolute mean value of the radial component. It corresponds in case of a synchronous generator to the temporal waveform and can therefore be compared to the formulation in (4). In this case, b and N equal the stator tooth pitch and number of teeth respectively. While in partial load condition, the fundamental oscillation of the tangential component is comparatively small, this finding cannot be verified in full load condition. The vibrational characterization with regard to magnetic force excitations has to consider both radial and tangential component in contrast to electrical machines with small air gap diameters [18]. This finding is confirmed by the results from the "Fast Fourier Transformation" (FFT) that are illustrated for partial and full load in Fig. 4. The temporal order l is normal- ized to the number of generator poles 2p. The difference between the amplitude of fundamental oscillation (order l / 2p = 1) of the radial and tangential between the radial and tangential force becomes significantly smaller in full load operation. The normalized orders 2 and 3 are moreover pronounced and might have impact to the overall vibrational characteristics that has to be studied by the force element. This model part is presented in the following chapter.

Multibody generator model
To analyze the generator, as one of the major vibration exciting element in the turbine, the generator model has been set up separately to identify critical vibrations of the generator. Structural components of the generator (stator, rotor, mounting structure, etc.) are modeled as elastic bodies by means of modal reduction as in the turbine model (see Fig. 5), whereas the nacelle cover, rotor blades and tower are left out. The electromagnetic forces are implemented via a novel force element in Simpack. Therefore, the air gap K surfaces of rotor and stator are divided into subsurfaces. Each subsurface has a marker in the multibody system, which are coupled to the elastic structure of the rotor and stator respectively via a kinematic coupling [13]. On each marker (rotor and stator) electromagnetic forces are exerted, that are derived from the magnetic flux density distribution in the air gap. The forces for each marker depend on the local air gap widths and the angular displacement between the stator marker and the rotor. These forces are pre-calculated by the means of a finite element (FE) model of the generator (see Sect. 3.1).
In Fig. 6a a schematic illustration of a pole pitch of the generator lamination is shown, with the angular displacement 1 between a stator marker (blue) and a rotor marker (yellow) and the local air gap width ı 1 . Fig. 6b also shows the corresponding radial forces, that are exerted onto the stator marker. During each simulation time step the local air gap width ı 1 and the angular displacement 1 is determined and the corresponding forces (radial and tangential) are exerted to the stator and rotor marker. This procedure is repeated for each marker and each simulation time step. By utilizing the symmetry of the generator, the precalculated FE solutions only have to include the forces for one pole pitch, instead of the whole circumference of the machine. Therefore, the amount of data needed from the FE simulation is decreased significantly.
The force characteristic for different angular displacements and air gap widths are simulated for several operating points along the power curve of the turbine. To model a continuous transition between the operating points, the electromagnetic forces are interpolated between the precalculated operating points, depending on the position on the power curve.

Results
With this model an exemplary test case can be created. In this test case the input torque of the drivetrain is ramped up constantly and the rotational speed changes according to the turbines power curve (see Figs. 7 and 8). The duration of the ramp-up is t ramp = 45s, with a sampling frequency of f sample 800 Hz. Included are 122 eigenmodes of the stator structure and 62 eigenmodes rotor structure, which are all modes up to f eigen;max 400 Hz. The runtime of this model has been approximately three hours.
The radial force distribution of one stator tooth (respectively stator marker) is shown in Fig. 9, which represent the force of one stator tooth. It can be seen, that the amplitude of the radial force is negative-which means it is directed into the direction of the air gap (closing the air gap). It can also be seen, that the stator tooth force fluctuates strongly, as the rotor angular position changes.
With the help of this ramp-up crucial eigenmodes can be identified, as the ramp-up covers the rotational-speed range of the generator. Figs. 10 and 11 show examples of the shape and acceleration of two different eigenmodes of the stator. The stator structure deforms in a shape of a hexagon or pentagon respectively.
It can clearly be seen, that in the ramp-up process the eigenmodes are excited during at particular rotational speeds. These rotational speeds correspond to the stator slot harmonic of the generator, which can be seen in Fig. 2. At t 30s the hexagonal shaped eigenmode is excited at a rotational speed of n 11 rpm. As the generator accelerates up to n 12.5 rpm the hexagonal oscillation decreases as the stator slot harmonic frequency drifts away from the eigenfrequency of f Hexagon . Subsequently the pentagonal eigenfrequency starts to oscillate as the stator slot harmonic drifts into the pentagonal eigenfrequency of f Pentagon .
Such a vibrational behavior may lead to acoustic emission of the turbine, depending on the acoustic transfer path from the generator to the emitting components (rotor blades, cover, etc.) of the turbine. However, with this model the generator behavior can be analyzed and improvements could be developed regarding the generator design or by adjusting the turbines behavior. E.g. the amount of stator slots could be changed or the wind turbine can be controlled in a way that these rotational speeds are not applied for stationary power production.

Conclusion and outlook
In this paper a multiphysical modelling approach for directdrive wind turbines is presented. This model is targeted towards simulating dynamic vibrations and structure-born acoustics of a gearless wind turbine application.
Due to the generator being one of the major excitation sources in this system, the generator is studied more thoroughly. The electromagnetic modelling of the generators flux density distribution is explained, which is fundamental for the magnetic forces on the rotor and stator surface. Furthermore, the multibody model of the drivetrain is introduced, in addition to an approach for applying Maxwell forces on the structural components of the generator.
The drivetrain behavior is studied for the situation of an exemplary ramp-up. During this case the input torque for the drivetrain is increased constantly and the generator air gap torque follows its torque-speed-characteristic. For this ramp-up the radial stator tooth force (marker force) is shown. These radial tooth forces excite the pentagonand hexagon-shaped eigenmodes of the generator structure, during the ramp-up.
It can be seen, that these eigenmodes are excited at particular rotational speeds, which correspond to the stator slot harmonic frequency.
This paper demonstrates, that this modelling approach can be used to simulated mechanical vibrations in gearless wind turbines, which helps to further optimize the vibrational and acoustic behavior of the turbines. In the now following research steps, the impact of wind loads onto the electromagnetic excitation will be examined, which is expected to enhance the electromagnetic excitation due to air gap widths change. Furthermore, a full turbine analysis coupled with the aerodynamic load simulation including rotor blades, tower, hub-and nacelle-cover will be conducted. For the electromagnetic modeling, future works will address the validation of the proposed 2D-modeling by 3Dsimulation. Particularly interesting aspects are the validation of the presented multi-slice approach as well as the study of 3D-air gap imperfections and their impact on axial force excitations. Moreover, additional force excitations by the rectifier feeding of the generator are under the scope of research.