FE-BE computation of electromagnetic noise of a permanent-magnetic excited synchronous ma-chine considering dynamic rotor eccentricity

Electromagnetic noise in Electrical Machines (EMs) occurs due to vibrations caused by magnetic forces acting onto rotor and stator surface. This is the dominant source for the considered permanent-magnetic excited synchronous machine in this paper. The radiated electromagnetic noise is sequentially calculated by a Finite Element (FE) and Boundary Element (BE) computation. An electromagnetic FE model is created to determine magnetic forces. Structure-borne sound and rotor dynamics are calculated using a structural dynamic FE model for the EM housing and the rotor. In order to predict resonance frequencies and amplitudes as reliable as possible, it is important to know the direction-dependent stiffness of the laminated rotor stacks and mechanical joints as well as their structural damping. Thereby, the properties of the laminated stack can be determined experimentally by a shear and dilatation test. Mechanical joint properties can be modelled by Thin-Layer Elements (TLEs) and the overall damping by the model of constant hysteretic damping. The radiated sound power is determined by a direct BE computation. The influence of dynamic rotor eccentricity on radiated sound power is examined for a run-up of the EM. All FE models are verified by data from experimental modal analysis.


Introduction and aim of this work
Acoustical noise of electrical machines is made up of three different types of sources: electromagnetically, mechanically and aerodynamically excited vibrations/noise [16].Electromagnetically caused vibrations result from electromagnetic forces in the air gap between rotor and stator [16].Mechanically caused vibrations are e.g.due to rotor dynamic loads, bearing defects or tolerances and shaft misalignments [16].Aerodynamic noise normally occurs due to the noise of cooling fans [16].In this paper, a Permanent Magnetic excited Synchronous Machine (PMSM) is considered with a totally enclosed water-cooled housing, see Fig. 1.The focus is set on electromagnetic noise.
For acoustical computations both, vibrations of the housing as well as of the rotor must be taken into account to model all structure-borne and airborne sound paths.Electromagnetic noise of electrical machines has already been considered in many publications.Most of them are only related to structure-borne sound / vibration analysis (FE for electromagnetics & structural dynamics) and sound radiation of the machine housing caused by electromagnetic force excitation on stator teeth (with/without eccentricity effects), cf.[10,22].Only in few works the rotor is included as additional component into the structural FE model (see [17,19,21]) and magnetic forces are applied on the rotor (see [19,21]).So far, no publication has been found by the authors which considers the numerical computation of rotor dynamics due to distributed electromagnetic forces on the laminated stack (from magnetic FE analysis) using three-dimensional (3D) structural solid FEs.Available works are based on analytical rotordynamic models with one or a few degrees of freedom or FE models built up by beam elements, cf.[1,21] and named references in [21].
Therefore, this work considers 3D FE rotordynamic analyses in order to examine electromagnetically caused rotor vibrations.The unknown direction-dependent stiffness and structural damping of the laminated stack are determined experimentally using a shear and dilatation test.In contrast to other approaches (named in [4]), this method allows for a physical implementation using only one experimental setup.For airborne sound predictions, a physical mechanical FE model of the housing is required for reliable predictions of structure-borne sound and vibration responses.In order to describe the damping of the system, the model of constant hysteretic damping can be efficiently applied, cf.[2].Local energy dissipation at mechanical interfaces is modelled by load-dependent parametrized TLEs [6,12] containing measured joint parameters, cf.[2,[13][14][15].To verify the FE models, Numerical Modal Analysis (NMA) results are compared to those of an Experimental Modal Analysis (EMA).The influence of rotor eccentricity on sound radiation is examined.

Parameter identification for laminated stacks and joints
In [4,5,7,8], a measurement set-up (Fig. 2) is presented for the determination of direction-dependent stiffness and structural damping of laminated rotor stacks which has been transferred and enhanced based on the original idea of Crandall et al. [9].The axial stiffness and structural damping can be determined by a dilatation test (Fig. 2(a)) [4,5].Similarly, the shear stiffness and structural damping can be derived from a shear test (Fig. 2(b)) [4,5].In both configurations, one test specimen is axially pre-tensioned between two adjacent steel plates by screw connections (chain lines in Fig. 2) according to the required prestress condition.The masses of the steel plates (including the distributed mass of the test specimen) are denoted by m1, m 2 and m 3 .Both set-ups are supported at the metal plates by two thin wires as indicated by crossed circular symbols.In order to examine the energy dissipation, both set-ups are excited harmonically by an electromagnetic shaker whereas the excitation direction is different by 90°.In a steady-state condition, the excitation force and acceleration of each metal plate are measured by a force sensor and accelerometers.The relative displacement Δx = x 3 -x 2 (dilatation test) and Δz i = z 1 -z i with i = 2, 3 (shear test) can be calculated by integrating the acceleration signals twice with respect to time and subtracting them from each other, cf.[2].Plotting the transmitted normal force  , =  2 • ̈2 versus the relative displacement Δx (dilatation test) or  ,, =   • ̈ versus Δz i with i = 2, 3 (shear test) leads to a hysteresis curve (cf.[15] and Fig. 2 The structural damping loss factors η x (dilatation test) and η z,i (shear test) can be computed using Eqs.( 1) by dividing the dissipative damping energy W D,x or W D,z,i (dilatation or shear test) through 2π times the maximum stored potential energy U max,x or U max,z,i , cf. [2,4,15].
For a certain pre-load F N , the normal stiffness can be derived from the slope of the hysteresis, cf.[15].Using the normal stiffness k x , initial thickness h 0 of the test specimen and its cross section A the Young's Modulus E can be computed by Eqs.(2).Similarly, the Shear Modulus G follows from shear stiffness k z,i , cf. [15].Further details can be found in [4,5].
In the 1980s, a generic lap joint resonator was developed [3] which allows for the determination of the tangential contact stiffness and structural damping of mechanical joints.This set-up has been enhanced in further projects [13,14,18]; also for gaskets.Details on parameter identification can be found in [2,3,6,[13][14][15].The joint damping loss factor strongly depends on the relative displacement, tangential force F T,x and normal force [6,12,14].This dependency must be taken into account in FE models by a load-dependent TLE parametrization for applications with inhomogeneous contact pressure distributions [6,12].In this work, joint parameters are taken from [13,14].

Numerical simulation
Based on the experimental identified parameters, a 3D ABAQUS TM FE model of the electrical machine housing and rotor is created.Electromagnetic forces computed with a 2D FE model (cf.also [1,16,17,21]) using the software FEMAG are distributed onto the stator teeth and the rotor surface and are Fourier analyzed according to [8] using MATLAB.The vibrations on the surface of the structure are computed by an individual frequency response analysis per electromagnetic order.The out-coming velocity field serves as Neumann Boundary Condition (BC) for the direct BE computation (exterior problem), cf.[20].
The FE model of the housing (cf.[12]) consists of all components shown in Fig. 1(b) while the rotor consists of all components shown Fig. 1(a).Mechanical interfaces and gaskets of the housing are modeled by TLEs using a load-dependent parametrization in order to take inhomogeneous contact pressure distributions into account [6,12].For the rotor, a homogeneous TLE parametrization is sufficient [7,8].Further details on FE models, contact modeling and contact pressure measurements can be found in [6,7,12].https://doi.org/10.1051/matecconf/201821118005VETOMAC XIV To compute eigenfrequencies, mode shapes and modal damping, a complex eigenvalue analysis of the housing and rotor model for a 'free-free' boundary condition is done.For harmonic vibrations in frequency domain, the complex eigenvalue problem is (cf.[12]) Thereby, M is the mass matrix, K * the complex stiffness matrix,  k the k-th complex eigenvalue and ϕ k the k-th complex eigenvector.Damping is considered by the model of constant hysteretic damping and the complex stiffness matrix given in Eq. ( 4) [2,12].Thereby, K is the real-valued stiffness matrix of the whole system, K c and η c are the realvalued stiffness matrix and loss factor of component c, K i and η i are the real-valued stiffness matrix and loss factor of the metallic interface i and K g and η g are the real-valued stiffness matrix and loss factor of gasket g.The joint stiffness K i of each interface is modeled by TLEs with an orthotropic elasticity matrix D i [2,[12][13][14][15].
The stiffness matrix K c of the laminated stack is modeled by a transversely isotropic material model [5] using parameters from the shear and dilation test [5,7,8].More details are given in [5,7].Values for η c (excepted the laminated stacks), η i and k T are taken from [13].Due to the asymmetric geometry, all rotordynamic analysis are carried out in a rotational coordinate system according to Eq. ( 5), cf.[11,23].C represents the Coriolis matrix, ′ the stiffness matrix enhanced by bearing stiffness, Z the centrifugal matrix,   the geometric stiffness matrix,  the displacement vector, ̇ the velocity vector, ̈ the acceleration vector and F(t) the excitation force vector (all in rotating coordinates!) [23].Ω is the angular rotor spin speed [23].
External viscous damping is neglected.Internal damping is modeled by Eq. ( 4) in frequency domain.Details on assumptions, damping and bearing models can be found in [7,8].In frequency domain, sound radiation can be described by the Helmholtz equation [20] ∆̂+  2 • ̂= 0 with  = /, where ̂ represents sound pressure,  circular wave number,  = 2 circular frequency, f frequency and c sound velocity [16,20,24].Considering the Sommerfeld radiation condition, describing boundary conditions (Dirichlet, Neumann or Rubin), method of weighted residuals combined with the collocation method and discretization with BEs onto the closed structural surface leads after certain mathematical operations to the frequency-dependent discretized boundary integral equation in matrix form [20,24]   () •  ̂ () =   () •  ̂  ()  () • () = () (7) in order to solve for the unknown boundary data () [20,24].In this work, the CHIEF (Combined Helmholtz Integral Equation Formulation) method is applied to reduce the influence of spurious modes [20,24].Based on computed boundary data, the radiated pressure field  ̂ () can be evaluated frequency-dependent using the collocation method applied to the representation formula leading to an algebraic equation system, Eq. ( 8) [20,24].Then, sound power is computed based on radiated sound pressure using ISO 3744 [25].

Computational results
Fig. 3(a) depicts a comparison of computed and measured modal damping factors using different TLE modeling approaches [6,12].It can be seen that the general damping behavior can be predicted quite well [12].Thereby, weakly damped modes (see Fig. 3(b), mode 4) due to slight energy dissipation at the mechanical interfaces as well as highly damped modes (Fig. 3(c)) can be identified correctly [12].Table 1 shows a comparison of a NMA and EMA for two equivalent built-up rotors with good agreement for the first few eigenfrequencies and mode shapes [7].Damping of the first two modes is predicted quite well but the error increases at higher frequencies [7,8].

Table 1.
Comparison of results of a numerical and experimental modal analysis [7,8].Fig. 4 depicts spectrograms of the computed radiated sound power for a run-up containing different time harmonic EM orders.It shows the effect of dynamical rotor eccentricity.In contrast to zero eccentricity, further orders (±1) arises, cf.[10,16].A resultant radial force (Unbalanced Magnetic Pull, UMP [11]) occurs additionally to torsional forces.The corresponding rotordynamic results for force order k=36 can be seen in Fig. 4, cf.[7,8].