Review of Electromagnetic Vibration in Electrical Machines

Electrical machines are important devices that convert electric energy into mechanical work and are widely used in industry and people’s life. Undesired vibrations are harmful to their safe operation. Reviews from the viewpoint of fault diagnosis have been conducted, while summaries from the perspective of dynamics is rare. This review provides systematic research outlines of this field, which can help a majority of scholars grasp the ongoing progress and conduct further investigations. This review mainly generalizes publications in the past decades about the dynamics and vibration of electrical machines. First the sources of electromagnetic vibration in electrical machines are presented, which include mechanical and electromagnetic factors. Different types of air gap eccentricity are introduced and modeled. The analytical methods and numerical methods for calculating the electromagnetic force are summarized and explained in detail. The exact subdomain analysis, magnetic equivalent circuit, Maxwell stress tensor, winding function approach, conformal mapping method, virtual work principle and finite element analysis are presented. The effects of magnetic saturation, slot and pole combination and load are discussed. Then typical characteristics of electromagnetic vibration are illustrated. Finally, the experimental studies are summarized and the authors give their thoughts about the research trends.


Introduction
With the continuous development of the economy, electrical machines have been widely used in industries and people's lives.Ever higher performance requirements are being put forward for electrical machines.The magnetic field can interact with mechanical structures, which will produce unbalanced magnetic forces and excite harmful vibrations.These forces may have significant effects on the dynamic behavior and noise of rotors.The vibration and noise range of electrical machines is one of the important indicators in the manufacture.National standards also clearly prescribe the vibration and noise limits of rotating electrical machines.Vibration and noise will affect people's daily life and severe vibration will cause significant economic losses.Furthermore, the vibrations of electrical machines may decrease the efficiency of the energy conversion because the vibration and possible related temperature rise are unwanted energy losses for the system.
The vibrations of electrical machines can be divided into three categories: mechanical vibrations, electromagnetic vibrations and aerodynamic vibrations.Benefitting from the continuous improvement of design and manufacturing level, the performance of electrical machines has been greatly improved and their volume has become very small.For the widely used small and medium-sized electrical machines, electromagnetic vibrations are the main type.With the growth of living standards, people will pay more and more attention to the vibration of electrical machines.Investigating the Energies 2018, 11, 1779 3 of 33 of electromagnetic vibration are basically summarized from the view of fault diagnosis [22][23][24][25][26][27], and summaries from the perspective of dynamics are rather insufficient.With the continuous progress of research, some new technologies and methods are emerging.It is necessary to generalize the latest research progress of electromagnetic vibration from the perspective of dynamics and vibration to avoid repetitive work.In addition, a review which provides systematic research outlines and references can be beneficial for the majority of scholars in this field to promote the ongoing progress and development of the investigations.
Different from the condition monitoring standpoint, this review is mainly concerned with the dynamic issues of electromagnetic vibrations in electrical machines.First, the vibration sources, which include the mechanical and electromagnetic aspects, are summarized in Section 2. Then Section 3 presents in detail different analytical and numerical calculation methods for the electromagnetic force.After that, the electromagnetic vibration characteristics and experimental investigations are demonstrated in Section 4. Finally Section 5 summarizes the authors' thoughts about the trends and prospects of this research.

Sources of Electromagnetic Vibration
Under ideal conditions, the air gap between the stator and rotor is uniform and the magnetic circuit is symmetrical.The rotor rotates in the uniform magnetic field and the total force of the radial electromagnetic force is zero.If mechanical or electromagnetic factors make the radial force around the rotor circumference uneven, it will produce an electromagnetic force which is also known as the unbalanced magnetic force (UMF).UMF will cause undesired electromagnetic vibration and noise, exacerbate the bearing wear, influence the stability of the rotor system and even produce rubbing between the rotor and stator.The potential hazards are prominent.Therefore, the investigation of this coupling interaction is necessary and important.
In essence, the main source of electromagnetic vibration is the asymmetry of the magnetic circuit in the electrical machine.The misalignment between the stator and rotor is the most common cause of UMF.Furthermore, the uneven magnetization of the material and the improper winding can also generate UMF.Generally speaking, the electromagnetic sources can be divided into mechanical factors and electromagnetic factors.

Mechanical Sources
The mechanical causes of the electromagnetic force are mainly the air gap eccentricity between the stator and the rotor.As Figure 1 indicates, the sources of air gap eccentricity can be divided into four categories: shape deviation, parallel eccentricity, inclined eccentricity and curved eccentricity.The surface corrugations of the outer rotor circle and inner stator circle will affect the uniformity of the air gap length.In addition, when the stator and rotor are not regular cylinders, their shape deviation can produce air gap eccentricity.Lundström et al. [28,29] studied the air gap eccentricity and electromagnetic force caused by the deviation of generator shape in detail.The characteristics of dynamic responses including the whirling frequency and amplitude were investigated.With the progress of mechanical manufacturing technology, the probability of shape deviation in the rotor and stator is gradually declining.Another very important and widely investigated air gap eccentricity type is that the rotor shaft does not coincide with the stator axis.It is assumed that the stator and rotor are ideal cylinders and have a smooth surface.Moreover the stator axis is straight.The air gap changes caused by assembly error and bearing wear etc. can be regarded as different air gap eccentricities.When the rotor shaft is straight, eccentricity can be further divided into parallel eccentricity and inclined eccentricity.When the rotor shaft is bending, the eccentricity type is curved eccentricity.

Parallel Eccentricity
As shown in Figure 2, when the rotor shaft and the stator axis are parallel, the air gap eccentricity can be divided into three categories: one is the static eccentricity which refers to that the air gap eccentricity already exists before operation and the rotor rotates with its own geometric center axis [30,31].The second category is the dynamic eccentricity which occurs when the stator and rotor are concentric at first and the eccentricity occurs during the operation.The rotor rotates with the geometric axis of the stator [32,33].Static eccentricity and dynamic eccentricity are the most basic eccentricity types.The third category is the static and dynamic mixed eccentricity, that is, static eccentricity and dynamic eccentricity coexist [33][34][35].The dynamic eccentricity is mainly caused by the mass unbalance of the rotor.The radial centrifugal force is generated during the rotation of the electrical machine, which results in the uneven air gap between the stator and the rotor.The static eccentricity is easily caused by the installation parallel deviation and bearing wear.Static and dynamic mixed eccentricity can be regarded as the static eccentricity plus the dynamic eccentricity.The air gap length for parallel eccentricity is of great importance to the calculation of the UMF and the vibration analysis of the rotor system.Taking the static and dynamic mixed eccentricity (as Figure 3 indicates) as an example, the air gap length formula is derived as follows [36]:

Parallel Eccentricity
As shown in Figure 2, when the rotor shaft and the stator axis are parallel, the air gap eccentricity can be divided into three categories: one is the static eccentricity which refers to that the air gap eccentricity already exists before operation and the rotor rotates with its own geometric center axis [30,31].The second category is the dynamic eccentricity which occurs when the stator and rotor are concentric at first and the eccentricity occurs during the operation.The rotor rotates with the geometric axis of the stator [32,33].Static eccentricity and dynamic eccentricity are the most basic eccentricity types.The third category is the static and dynamic mixed eccentricity, that is, static eccentricity and dynamic eccentricity coexist [33][34][35].The dynamic eccentricity is mainly caused by the mass unbalance of the rotor.The radial centrifugal force is generated during the rotation of the electrical machine, which results in the uneven air gap between the stator and the rotor.The static eccentricity is easily caused by the installation parallel deviation and bearing wear.Static and dynamic mixed eccentricity can be regarded as the static eccentricity plus the dynamic eccentricity.

Parallel Eccentricity
As shown in Figure 2, when the rotor shaft and the stator axis are parallel, the air gap eccentricity can be divided into three categories: one is the static eccentricity which refers to that the air gap eccentricity already exists before operation and the rotor rotates with its own geometric center axis [30,31].The second category is the dynamic eccentricity which occurs when the stator and rotor are concentric at first and the eccentricity occurs during the operation.The rotor rotates with the geometric axis of the stator [32,33].Static eccentricity and dynamic eccentricity are the most basic eccentricity types.The third category is the static and dynamic mixed eccentricity, that is, static eccentricity and dynamic eccentricity coexist [33][34][35].The dynamic eccentricity is mainly caused by the mass unbalance of the rotor.The radial centrifugal force is generated during the rotation of the electrical machine, which results in the uneven air gap between the stator and the rotor.The static eccentricity is easily caused by the installation parallel deviation and bearing wear.Static and dynamic mixed eccentricity can be regarded as the static eccentricity plus the dynamic eccentricity.The air gap length for parallel eccentricity is of great importance to the calculation of the UMF and the vibration analysis of the rotor system.Taking the static and dynamic mixed eccentricity (as Figure 3 indicates) as an example, the air gap length formula is derived as follows [36]: The air gap length for parallel eccentricity is of great importance to the calculation of the UMF and the vibration analysis of the rotor system.Taking the static and dynamic mixed eccentricity (as Figure 3 indicates) as an example, the air gap length formula is derived as follows [36]: A two-dimensional Cartesian orthogonal coordinate system is established.Os is the geometric center of the stator.O is initial geometric center of the rotor and Or is geometric center of the rotor during operation.By geometric relationship derivation, the air-gap length for static and dynamic mixed eccentricity can be approximately expressed as: where δ0 is the average air-gap length when the rotor is centered.r0 is the static eccentricity and r is the dynamic eccentricity.α is the air-gap angle with respect to x-axis.φ and θ are the angles of static eccentricity and dynamic eccentricity with reference to x-axis, respectively.If just the static eccentricity or dynamic eccentricity exists, Equation (1) can be respectively simplified as:

Inclined Eccentricity
In engineering practice, the height difference between the bearings on both sides of the rotor shaft and the inclination of the shaft etc. will cause air gap differences in the rotor at different positions along the axial direction.This eccentricity is named inclined eccentricity.In 1992 Akiyama [37] proposed inclined eccentricity based on the actual engineering needs.The inclined eccentricity can be further subdivided into symmetrical inclined eccentricity and mixed inclined eccentricity.Symmetric inclined eccentricity means that there is only angular deviation between the stator axis and the rotor shaft.For the mixed inclined eccentricity, there exist angular deviation and radial displacement between the stator axis and the rotor shaft.The mixed inclined eccentricity can be considered as a combination of symmetric inclined eccentricity and parallel eccentricity.Figure 4 shows the schematic diagram of symmetrical inclined eccentricity.Due to the fact that the A two-dimensional Cartesian orthogonal coordinate system is established.O s is the geometric center of the stator.O is initial geometric center of the rotor and O r is geometric center of the rotor during operation.By geometric relationship derivation, the air-gap length for static and dynamic mixed eccentricity can be approximately expressed as: where δ 0 is the average air-gap length when the rotor is centered.r 0 is the static eccentricity and r is the dynamic eccentricity.α is the air-gap angle with respect to x-axis.ϕ and θ are the angles of static eccentricity and dynamic eccentricity with reference to x-axis, respectively.If just the static eccentricity or dynamic eccentricity exists, Equation (1) can be respectively simplified as:

Inclined Eccentricity
In engineering practice, the height difference between the bearings on both sides of the rotor shaft and the inclination of the shaft etc. will cause air gap differences in the rotor at different positions along the axial direction.This eccentricity is named inclined eccentricity.In 1992 Akiyama [37] proposed inclined eccentricity based on the actual engineering needs.The inclined eccentricity can be further subdivided into symmetrical inclined eccentricity and mixed inclined eccentricity.Symmetric inclined eccentricity means that there is only angular deviation between the stator axis and the rotor shaft.For the mixed inclined eccentricity, there exist angular deviation and radial displacement between the stator axis and the rotor shaft.The mixed inclined eccentricity can be considered as a combination of symmetric inclined eccentricity and parallel eccentricity.Figure 4 shows the schematic diagram of symmetrical inclined eccentricity.Due to the fact that the air gap length of inclined eccentricity is related to the inclined angle and axial position, it is necessary to investigate this problem in the three-dimensional coordinate system.The air gap length of mixed inclined eccentricity is further investigated in the following paragraphs.The general case of air gap changes along the axial direction can be considered as mixed inclined eccentricity [38,39].As Figure 5 demonstrates, it is assumed that the stator is stationary and the midpoint of rotor in the axial direction is selected to be the origin of the coordinates.There exist two orthogonal coordinate systems in the three-dimensional space and they are parallel.One is the stationary coordinate system O-XYZ for the stator.The other is the unfixed coordinate system o-xyz for the rotor [40].The parameters that describe the rotor state can be divided into two categories: one is the dynamic parameters (x, y and θx, θy) describing the dynamic displacement and angle responses respectively, and the other is static parameters (x0, y0 and θx0, θy0) representing the static displacement and angle eccentricities respectively.The dynamic parameters are determined by the dynamic responses of the rotor system.However, the static parameters need to be provided in the beginning.
The complex three-dimensional problem can be decomposed to be a two-dimensional eccentricity by cutting the axis into numerous cross sections along the axial direction.The comprehensive eccentricity of the investigated cross section is: The air gap length of mixed inclined eccentricity is further investigated in the following paragraphs.The general case of air gap changes along the axial direction can be considered as mixed inclined eccentricity [38,39].As Figure 5 demonstrates, it is assumed that the stator is stationary and the midpoint of rotor in the axial direction is selected to be the origin of the coordinates.There exist two orthogonal coordinate systems in the three-dimensional space and they are parallel.One is the stationary coordinate system O-XYZ for the stator.The other is the unfixed coordinate system o-xyz for the rotor [40].
Energies 2018, 11, x 6 of 34 air gap length of inclined eccentricity is related to the inclined angle and axial position, it is necessary to investigate this problem in the three-dimensional coordinate system.The air gap length of mixed inclined eccentricity is further investigated in the following paragraphs.The general case of air gap changes along the axial direction can be considered as mixed inclined eccentricity [38,39].As Figure 5 demonstrates, it is assumed that the stator is stationary and the midpoint of rotor in the axial direction is selected to be the origin of the coordinates.There exist two orthogonal coordinate systems in the three-dimensional space and they are parallel.One is the stationary coordinate system O-XYZ for the stator.The other is the unfixed coordinate system o-xyz for the rotor [40].The parameters that describe the rotor state can be divided into two categories: one is the dynamic parameters (x, y and θx, θy) describing the dynamic displacement and angle responses respectively, and the other is static parameters (x0, y0 and θx0, θy0) representing the static displacement and angle eccentricities respectively.The dynamic parameters are determined by the dynamic responses of the rotor system.However, the static parameters need to be provided in the beginning.

Stator
The complex three-dimensional problem can be decomposed to be a two-dimensional eccentricity by cutting the axis into numerous cross sections along the axial direction.The comprehensive eccentricity of the investigated cross section is: The parameters that describe the rotor state can be divided into two categories: one is the dynamic parameters (x, y and θ x , θ y ) describing the dynamic displacement and angle responses respectively, and the other is static parameters (x 0 , y 0 and θ x0 , θ y0 ) representing the static displacement and angle eccentricities respectively.The dynamic parameters are determined by the dynamic responses of the rotor system.However, the static parameters need to be provided in the beginning.
The complex three-dimensional problem can be decomposed to be a two-dimensional eccentricity by cutting the axis into numerous cross sections along the axial direction.The comprehensive eccentricity of the investigated cross section is: where the value of the coordinate Z is the intersection point between the investigated cross section and OZ axis.
Energies 2018, 11, 1779 The common range of ϕ is [0, 2π], while the scope of inverse trigonometric functions is [0, π].Hence, the following extension is conducted: The air-gap length is a function of the air-gap angle, time and axial position [41][42][43].The unified air-gap length in an arbitrary cross-section can be approximately expressed by the equation as follows: The effective air-gap length along the axial direction becomes short because the rotor is inclined with respect to the stator.The actual air-gap interaction length between the rotor and the stator is determined by the static angle eccentricities and the dynamic angle responses: where L is the axial length of the air-gap.

Curved Eccentricity
The curved eccentricity always occurs to some extent in most large motors where the axis is bending.For example, a three-phase diving induction motor holds dynamic arc eccentricity.Due to the effects of load or the insufficient shaft stiffness, the shaft will also bend and form curved air gap eccentricity [44][45][46].The usual way to analyze this complex situation is to treat the electrical rotor as a number of small slices.In each slice, the air gap length can be analyzed according to the pattern of the basic eccentricities (parallel eccentricity or inclined eccentricity).
As shown in Figure 6, in order to obtain the air gap length at different positions, a multi-layer model is designed.The rotor is divided into many layers along the axial direction.Each layer is small enough so that the parallel eccentricity or inclined eccentricity can be applied [47].
where the value of the coordinate Z is the intersection point between the investigated cross section and OZ axis.
The common range of φ is [0, 2π], while the scope of inverse trigonometric functions is [0, π].Hence, the following extension is conducted: The air-gap length is a function of the air-gap angle, time and axial position [41][42][43].The unified air-gap length in an arbitrary cross-section can be approximately expressed by the equation as follows: The effective air-gap length along the axial direction becomes short because the rotor is inclined with respect to the stator.The actual air-gap interaction length between the rotor and the stator is determined by the static angle eccentricities and the dynamic angle responses: where L is the axial length of the air-gap.

Curved Eccentricity
The curved eccentricity always occurs to some extent in most large motors where the axis is bending.For example, a three-phase diving induction motor holds dynamic arc eccentricity.Due to the effects of load or the insufficient shaft stiffness, the shaft will also bend and form curved air gap eccentricity [44][45][46].The usual way to analyze this complex situation is to treat the electrical rotor as a number of small slices.In each slice, the air gap length can be analyzed according to the pattern of the basic eccentricities (parallel eccentricity or inclined eccentricity).
As shown in Figure 6, in order to obtain the air gap length at different positions, a multi-layer model is designed.The rotor is divided into many layers along the axial direction.Each layer is small enough so that the parallel eccentricity or inclined eccentricity can be applied [47].To sum up, the common types of air gap eccentricity and logical relationship between them can be summarized in Figure 7:

Stator
To sum up, the common types of air gap eccentricity and logical relationship between them can be summarized in Figure 7

Electromagnetic Sources
The sources in electromagnetic aspects can be summarized in four categories: short circuit, open circuit, magnetization unevenness and winding topology asymmetry.Under normal circumstances, the air gap flux and the electromagnetic force distribution are even and symmetrical.When a short circuit occurs in the rotor or stator slot, the magnetic flux of air gap changes.The UMF results in radial vibration [48,49].The reasons for short circuits include the unfixed excitation winding end, winding deformation, winding manufacturing defects and foreign matter intake, etc.The current flowing through the short-circuited coil is zero, which causes a decrease in the magnetic potential of the corresponding magnetic pole and an asymmetry of the magnetic field.Thereby the UMF is generated.The common types of short circuit are turn-to-turn, coil-to coil, phase-to-phase and phase-to ground short-circuit [23].Wan et al. [50] studied the influence of the short circuit on the force acting on the rotor and pointed out that the inter-turn short circuit in generator rotor caused thermal imbalance and magnetic imbalance.Wallin et al. [51] studied the UMF and flux density distribution resulting from winding inter turn-short circuits.Through experimental and numerical simulations, it was found that an additional unbalanced electromagnetism was generated.
The open circuit can also make the magnetic circuit unequal [52,53].Broken rotor bars and broken end rings are the most common open circuit types in a squirrel-cage rotor [54,55].If the broken bars are distributed over the poles, the current of the broken bars flows to the adjacent bars.This leads to unbalanced magnetic flux.If the broken bars are adjacent to each other, the current of each broken bar may not flow to its adjacent bar, therefore a more uneven magnetic flux distribution occurs.Therefore both the adjacent broken bars and distributed broken bars can result in UMF [56].Based on the electromagnetic theory, Jung et al. [57] derived a corrosion model for a rotor-bar-fault induction motors.In addition, Baccarini et al. [58] proposed an analytical model for induction machines considered the broken rotor bars and other factors.
Moreover, when the magnetization of the motor material is not uniform, the electromagnetic force per unit area at the rotor surface is different.Therefore the electromagnetic force is not zero [59].The factors that cause uneven magnetization can be divided into non-uniform magnetization of permanent magnet materials and magnetization of soft magnetic materials.The uniformity of the magnetization is affected by many factors, such as the aging of the magnet, the mutual repulsion of the magnetic field, the asymmetric magnetization during manufacture and the magnetic edge effect of the magnetic ring etc.The magnetization inhomogeneity will make the magnetic circuit asymmetric and produce the electromagnetic force.In addition, the asymmetry winding of the electrical machine will also cause electromagnetic force [60][61][62][63].Zhu et al. [64]

Electromagnetic Sources
The sources in electromagnetic aspects can be summarized in four categories: short circuit, open circuit, magnetization unevenness and winding topology asymmetry.Under normal circumstances, the air gap flux and the electromagnetic force distribution are even and symmetrical.When a short circuit occurs in the rotor or stator slot, the magnetic flux of air gap changes.The UMF results in radial vibration [48,49].The reasons for short circuits include the unfixed excitation winding end, winding deformation, winding manufacturing defects and foreign matter intake, etc.The current flowing through the short-circuited coil is zero, which causes a decrease in the magnetic potential of the corresponding magnetic pole and an asymmetry of the magnetic field.Thereby the UMF is generated.The common types of short circuit are turn-to-turn, coil-to coil, phase-to-phase and phase-to ground short-circuit [23].Wan et al. [50] studied the influence of the short circuit on the force acting on the rotor and pointed out that the inter-turn short circuit in generator rotor caused thermal imbalance and magnetic imbalance.Wallin et al. [51] studied the UMF and flux density distribution resulting from winding inter turn-short circuits.Through experimental and numerical simulations, it was found that an additional unbalanced electromagnetism was generated.
The open circuit can also make the magnetic circuit unequal [52,53].Broken rotor bars and broken end rings are the most common open circuit types in a squirrel-cage rotor [54,55].If the broken bars are distributed over the poles, the current of the broken bars flows to the adjacent bars.This leads to unbalanced magnetic flux.If the broken bars are adjacent to each other, the current of each broken bar may not flow to its adjacent bar, therefore a more uneven magnetic flux distribution occurs.Therefore both the adjacent broken bars and distributed broken bars can result in UMF [56].Based on the electromagnetic theory, Jung et al. [57] derived a corrosion model for a rotor-bar-fault induction motors.In addition, Baccarini et al. [58] proposed an analytical model for induction machines considered the broken rotor bars and other factors.
Moreover, when the magnetization of the motor material is not uniform, the electromagnetic force per unit area at the rotor surface is different.Therefore the electromagnetic force is not zero [59].The factors that cause uneven magnetization can be divided into non-uniform magnetization of permanent magnet materials and magnetization of soft magnetic materials.The uniformity of the magnetization is affected by many factors, such as the aging of the magnet, the mutual repulsion of the magnetic field, the asymmetric magnetization during manufacture and the magnetic edge effect of the magnetic ring etc.The magnetization inhomogeneity will make the magnetic circuit asymmetric and produce the electromagnetic force.In addition, the asymmetry winding of the electrical machine Energies 2018, 11, 1779 9 of 33 will also cause electromagnetic force [60][61][62][63].Zhu et al. [64] analyzed the electromagnetic force characteristics for permanent magnet motor in the case of completely asymmetrical winding.

Calculation Method of the Electromagnetic Force
The calculation of UMF is an important part of any electromagnetic vibration analysis.For the sake of dynamic modeling, the UMF is usually decomposed into a radial force directed to the shortest air gap and a tangential force directed perpendicular to the radial one.Based on the assumption that the magnitude of the electromagnetic force is proportional to eccentricities, Tenhunen et al. [65][66][67] studied the radial and tangential electromagnetic forces of the eccentric rotor at different rotational speeds.Frosini et al. [68] investigated the effects of eccentricity on the radial and tangential electromagnetic forces at open circuit and load.The authors established an analytical function of the electromagnetic force with respect to the known parameters.Wu et al. [61] studied the radial and tangential electromagnetic forces of surface-mounted permanent magnet motors under load and revealed the mechanism of increasing and decreasing the radial and tangential forces.Since the radial component of the air gap flux is much larger than the tangential component, the tangential component is generally neglected in the calculation [69][70][71][72][73][74].
As Figure 8 displays, the calculation method of UMF has undergone a complex development process.The early research mainly focused on the theoretical analysis of electromagnetic force and the linear expression was adopted [75][76][77].Behrend [75] obtained a linear formula of the electromagnetic force based on the hypothesis that the UMF is proportional to the air gap eccentricity.Calleecharan et al. [76] simplified the electromagnetic force as a linear spring which holds a negative stiffness coefficient when studying industrial hydro-generators.Although linear expressions are simple and convenient, the results are only reliable when the air gap eccentricity is small enough.In 1965, Funke et al. [11] suggested that there was a nonlinear function between the electromagnetic force and the air gap eccentricity.The nonlinear calculation methods of the UMF were widely studied since then.
Energies 2018, 11, x 9 of 34 analyzed the electromagnetic force characteristics for permanent magnet motor in the case of completely asymmetrical winding.

Calculation Method of the Electromagnetic Force
The calculation of UMF is an important part of any electromagnetic vibration analysis.For the sake of dynamic modeling, the UMF is usually decomposed into a radial force directed to the shortest air gap and a tangential force directed perpendicular to the radial one.Based on the assumption that the magnitude of the electromagnetic force is proportional to eccentricities, Tenhunen et al. [65][66][67] studied the radial and tangential electromagnetic forces of the eccentric rotor at different rotational speeds.Frosini et al. [68] investigated the effects of eccentricity on the radial and tangential electromagnetic forces at open circuit and load.The authors established an analytical function of the electromagnetic force with respect to the known parameters.Wu et al. [61] studied the radial and tangential electromagnetic forces of surface-mounted permanent magnet motors under load and revealed the mechanism of increasing and decreasing the radial and tangential forces.Since the radial component of the air gap flux is much larger than the tangential component, the tangential component is generally neglected in the calculation [69][70][71][72][73][74].
As Figure 8 displays, the calculation method of UMF has undergone a complex development process.The early research mainly focused on the theoretical analysis of electromagnetic force and the linear expression was adopted [75][76][77].Behrend [75] obtained a linear formula of the electromagnetic force based on the hypothesis that the UMF is proportional to the air gap eccentricity.Calleecharan et al. [76] simplified the electromagnetic force as a linear spring which holds a negative stiffness coefficient when studying industrial hydro-generators.Although linear expressions are simple and convenient, the results are only reliable when the air gap eccentricity is small enough.In 1965, Funke et al. [11] suggested that there was a nonlinear function between the electromagnetic force and the air gap eccentricity.The nonlinear calculation methods of the UMF were widely studied since then.In recent years, with the increasing reliability requirements for motor product, the nonlinear calculation of electromagnetic force has attracted extensive attention.Many calculation methods are In recent years, with the increasing reliability requirements for motor product, the nonlinear calculation of electromagnetic force has attracted extensive attention.Many calculation methods are Energies 2018, 11, 1779 10 of 33 emerging and generally there are three major categories: analytic methods, numerical methods and combinations of analytical and numerical methods.

Exact Subdomain Analysis
Exact subdomain analysis (ESA) is a method that divides the whole domain of the electrical machine into several subdomains and each subdomain is solved precisely.As Figure 9 displays, the solution domain in the ESA method is divided into five parts: air gap domain, stator core domain, rotor core domain, outer boundary domain and permanent magnet domain.There are boundary conditions between the regions and the radial and tangential components of the electromagnetic force can be considered.To obtain the solution of the air-gap field distribution, some basic assumptions for ESA are needed.For instance, the permeability of stator/rotor iron is infinite and the saturation as well as end effect are neglected.
Energies 2018, 11, x 10 of 34 emerging and generally there are three major categories: analytic methods, numerical methods and combinations of analytical and numerical methods.

Exact Subdomain Analysis
Exact subdomain analysis (ESA) is a method that divides the whole domain of the electrical machine into several subdomains and each subdomain is solved precisely.As Figure 9 displays, the solution domain in the ESA method is divided into five parts: air gap domain, stator core domain, rotor core domain, outer boundary domain and permanent magnet domain.There are boundary conditions between the regions and the radial and tangential components of the electromagnetic force can be considered.To obtain the solution of the air-gap field distribution, some basic assumptions for ESA are needed.For instance, the permeability of stator/rotor iron is infinite and the saturation as well as end effect are neglected.The scalar magnetic potential distribution in the air gap, stator iron, rotor iron and the exterior region is governed by the Laplace equation.While the scalar magnetic potential distribution in the magnet domain is governed by the quasi-Poissonian equation.These equations are expressed as follows: The scalar magnetic potential distribution in the air gap, stator iron, rotor iron and the exterior region is governed by the Laplace equation.While the scalar magnetic potential distribution in the magnet domain is governed by the quasi-Poissonian equation.These equations are expressed as follows: where r and θ are the radius and angle of the investigated subdomain in the polar coordinates, respectively.ϕ A , ϕ M , ϕ S , ϕ R and ϕ O represent the magnetic scalar potential in the air gap, magnet, stator, rotor and the exterior region, respectively.Equations ( 8)-( 12) are usually solved by variable separation technique and many scholars have applied the SEA method to investigate the magnetic field distribution in the electrical machine.The ESA method works for radial or parallel magnetized magnets as well as for the overlapping or nonoverlapping stator windings [78].In addition, the slotting effects [79,80] and different slot and pole combinations [81] can be considered in the SEA model.In a specific research, the five subdomains are often simplified.For instance, two domains including the air gap domain and permanent magnet domain were investigated [82].Three domains consist of magnet, air gap and slots were developed [83,84].Moreover, the ESA method is applicable to many types of electrical machines such as brushless permanent magnet machines and surface mounted permanent-magnet machines etc. Rahideh et al. [14] established a polar coordinate system to analyze the air gap/winding field and permanent magnet field distribution of slotless brushless permanent magnet motor with open circuit.Kumar et al. [85] proposed an improved analytical model by developing the instantaneous air-gap field distribution for a permanent magnet brushless DC motor.

Magnetic Equivalent Circuit
As Figure 10 illustrates, in the magnetic equivalent circuit (MEC) method, the rotor yoke, rotor teeth, rotor leakage, air gap, stator leakage, stator teeth and stator yoke are equivalent to a voltage loop and are superimposed along the magnetic circuit.The nodes of these voltage loop represent scalar magnetic potentials at different positions in the electrical machines.The current in the voltage loop passes through each node and it means magnetic flux goes through the magnetic unit.Figure 10 is an illustration of a permanent-magnet synchronous machine.The physical connection between the left and right ends of the magnetic circuit is modeled through common variables.A linear leakage permeance exists between every two stator teeth due to stator slot opening.The magnetic equivalent circuit method is mainly based on Kirchhoff's law and Gauss's law.The equations are solved by the Gaussian elimination method.The MEC can be regarded as a compromise between the electrical lumped-parameter models and finite element analysis.This approach shows the advantage of a close relationship with the physical field distributions in electrical machines.Moderate complexity and reasonable accuracy are reflected in the calculation.When compared with the finite element method, the disadvantages may be that the eddy current and skin effect cannot be handled perfectly.
The MEC method was formally introduced to electromechanical systems by Laithwaite [86] and Carpenter [87], respectively.In the late 1980s, Ostovic developed a series of publications about MEC modeling of induction and permanent-magnet synchronous machines [88][89][90][91].Since then, the application of MEC method was extensively investigated.The publications are mainly about electrical machines and selected examples are as follow: Sudhoff et al. [92] used the MEC method for induction motor modeling and proposed a methodology of constructing a state-variable model.Serri et al. [93] applied the MEC method to analyze the torque and radial force of a multi-phase bearingless permanent magnet motor, which simplified the design process of the motor.Xie et al. [94] studied the air-gap flux density of a dual-rotor permanent magnet motor by the MEC method and discussed the effects of winding current harmonics, winding types and radial and parallel magnetization types.Based on Kirchhoff's second law of magnetic circuit, Hanic et al. [95] proposed an analysis approach for saturated surface-mounted permanent magnet motor with no-load by using the conformal mapping and MEC method.Fleming et al. [96] conducted real-time simulation of switched reluctance motor by MEC.
mapping and MEC method.Fleming et al. [96] conducted real-time simulation of switched reluctance motor by MEC.

Maxwell Stress Tensor
The Maxwell stress tensor (MST) method describes the interaction between current and magnetic field as well as the distribution of magnetic flux density over the contour of the body.The basic laws of macroscopic electromagnetic phenomena can be expressed by the Maxwell equations.The total force on an object can be obtained by integrating the Maxwell stress tensor over a closed surface enclosing the object [98].Therefore, when the air gap in electrical machine is considered, the integration surface can be the circular plane between the two bodies and has a cylindrical shape with the normal vector pointing outwards.The Maxwell stress tensor  is utilized to calculate electromagnetic force F on a moving body: where volume V contains the investigated object.The Maxwell stress tensor is independent of coordinate system and defined as: where B is the magnetic flux density.i, j represent the components in the specific coordinate system and δij is the Kronecker delta function.By considering the Gauss's theorem, the Maxwell stress tensor force can be rewritten in a more convenient form:

Maxwell Stress Tensor
The Maxwell stress tensor (MST) method describes the interaction between current and magnetic field as well as the distribution of magnetic flux density over the contour of the body.The basic laws of macroscopic electromagnetic phenomena can be expressed by the Maxwell equations.The total force on an object can be obtained by integrating the Maxwell stress tensor over a closed surface enclosing the object [98].Therefore, when the air gap in electrical machine is considered, the integration surface can be the circular plane between the two bodies and has a cylindrical shape with the normal vector pointing outwards.The Maxwell stress tensor K is utilized to calculate electromagnetic force F on a moving body: where volume V contains the investigated object.The Maxwell stress tensor is independent of coordinate system and defined as: where B is the magnetic flux density.i, j represent the components in the specific coordinate system and δ ij is the Kronecker delta function.By considering the Gauss's theorem, the Maxwell stress tensor force can be rewritten in a more convenient form: where S is a surface enclosing the investigated body and n is a unit vector normal to the surface.For the specific application of MST method in electrical machines, the basic idea is to obtain the Maxwell stress on the rotor surface and the detailed process is presented in the following.
Energies 2018, 11, 1779 13 of 33 The air-gap permeance can be calculated as: where µ 0 is the vacuum permeability The magnetic flux density distribution of the air-gap is [99,100]: where F(α,t) is the synthesis fundamental magnetomotive force (MMF).
The normal component and tangential component of Maxwell stress are: In general, the tangential component of the flux density is much smaller than the normal component and can be considered negligible.Then the Maxwell stress perpendicular to the core/air boundary is given by: The two common approaches calculating the Maxwell stress are the Fourier series method [101] and direct integral method [36].In the Fourier series method, the air-gap permeance can be expanded as a Fourier series: where ε is the relative eccentricity.
It can be figured out that the first three harmonic components are dominant for ordinary eccentricities.Therefore, by ignoring the higher permeance harmonics and making some simplifications, the resulting electromagnetic forces in the horizontal and vertical direction are obtained as follows: where p is the number of pole-pair and ω is the supply frequency of the electrical machine.It can be observed that the results change with polyphase excitations: R and L are the radius and length of the rotor, respectively.For the direct integration method, the UMF can be obtained by integration of the Maxwell stress on the rotor surface.The expressions for parallel eccentricity are: If the inclined eccentricity is taken into consideration, not only the electromagnetic force but also the electromagnetic torque need to be calculated: where l is the axial length of the air-gap and l satisfies the following equation: The last five years have witnessed some progress in the MST method.Meessen et al. [102] selected the MST method to calculate the magnetic force components in the cylindrical coordinate system.By inserting the analytical expressions, the method can be fast and accurate.Spargo et al. [103] developed a semi-numerical method to calculate the harmonic torque components based on the MST theory which provides a simple algebraic expression.Bermúdez et al. [104] extended the MST method to consider the nonlinear magnetic media and local force distribution.The resultant electromagnetic force was verified well.

Winding Function Approach
The winding function approach (WFA) takes all harmonics of the magnetomotive force into consideration, with no restrictions concerning the symmetry of stator windings and rotor bars.However, the classic WFA is not suitable for the modeling of eccentricities since it cannot consider air-gap variations, although it was initially also applicable of these cases.The mutual inductances between stator phase and rotor loops (L sr ) are different from those between the rotor loops and stator phase (L rs ), and it is difficult to find a physical meaning of this asymmetry.As a result, the modified winding function approach (MWFA) for the inductance calculation considering air-gap eccentricity was proposed [105].In this method, the air-gap constant is replaced by an air-gap function which depends on the relative position of the rotor with respect to the stator.This method has been applied to analyze static, dynamic and mixed eccentricity in induction machines [106,107].In addition, this modification can be further extended to consider axial skewing.
The basic idea of WFA can be explained in an induction machine.Taking an induction machine which has m stator circuits and n rotor bars as an example, the cage can be regarded as n identical and equally spaced rotor loops.Voltage equations for the induction machine can be written in vector-matrix form as follows [108,109]: where: The stator and rotor flux linkages are given by: L ss is an m × m matrix with the stator self and mutual inductances, L rr is an n × n matrix with the rotor self and mutual inductances, L sr is an m × n matrix composed by the mutual inductance between the stator phases and the rotor loops.L rs is an n × m matrix and L rs = L sr T .The mechanical equations for electrical machines are obtained: where θ r is the rotor position, ω is the angular speed and J rl is the rotor-load inertia.T e is the electromagnetic torque and T l is the load torque.The magnetic co-energy which stores in the magnetic circuits and can be written as: The electromagnetic torque T e can be obtained from the magnetic co-energy: According to the WFA, the general expression of the mutual inductance between any winding i and j in any electrical machine is given by: where g −1 (ϕ, θ) is the inverse air-gap length.N i (ϕ, θ) and N j (ϕ, θ) are the winding functions of the windings i and j. ϕ is the angle along the inner surface of the stator and θ is the angular position of the rotor with respect to the stator.The winding function N i (ϕ, θ) is a function of ϕ and θ for a rotating coil.While it is only a function of ϕ for a stationary coil.This expression is inappropriate to handle arbitrary distribution windings of synchronous machines and the MWFA should be employed.
Much work has been done on the WFA method.Faiz et al. extended the winding function theory for non-uniform air gap eccentricity in rotating electric machinery [110] and applied the MWFA to calculate the time-dependent inductances of the motor with static, dynamic and mixed eccentricities in a unified manner [111].Ghoggal et al. proposed a MWFA by taking account the skewing rotor bars effects [112] and teeth saturation due to local air-gap flux concentration [113].Tu et al. [114] investigated the actual winding arrangement of a synchronous machine based on the WFA and Iribarnegaray et al. [115] gave a critical review of the modified winding function theory.

Conformal Mapping Method
The Laplace equations of the magnetic field distribution in the electrical machine are hard to solve directly.A possible solution is to convert this problem into an orthogonal coordinate system by the conformal mapping method (CMM).The CMM can maintain the solution of Laplace's equation in both the original and transformed domain.If the field distribution (e.g., rectangle, circle) of a geometrical subdomain is known, the field distribution in more complex geometrical subdomains (e.g., slotted air gap) can be calculated by the CMM.The permeance of the electrical machine was often obtained by applying the unit magnetic potential in the CMM.The magnet is ignored and infinitely deep rectilinear stator slots are assumed.The results can be acceptable when the ratio of the air gap length to slot pitch is relatively high.Such a single slot model is appropriate for the electrical machine with a relatively small number of slots.However an electrical machine with plenty of slots are common in ordinary cases.Considering that the effect of slotting on the variation of air gap is similar in an electrical machine, the single-slot approach is extensively investigated [116].
As illustrated in Figure 11, the conformal transformation is applied to transform the geometric shape into a slotless air gap in which the field distribution can be solved.The solution is then mapped back into the complex plane where the actual slot shape exists.Four conformal transformations are required to transform the slotted air gap into a slotless air gap: The original geometry represents a single slot opening in the S plane.This geometric shape needs to be transformed into its linear model in the Z plane utilizing a logarithmic conformal transformation defined as [117,118]: where j s m jn re     , z x jy   .The relationship between the coordinate in the S and Z planes is: The next step is to convert the geometric structure in the Z plane into the upper half of the W plane through Schwarz-Christoffel transformation.The approaches of W plane to T plane and T plane to K plane are similar to that from Z plane to W plane.The individual transformations between the planes can be written as follows: The original geometry represents a single slot opening in the S plane.This geometric shape needs to be transformed into its linear model in the Z plane utilizing a logarithmic conformal transformation defined as [117,118]: where s = m + jn = re jθ , z = x + jy.The relationship between the coordinate in the S and Z planes is: The next step is to convert the geometric structure in the Z plane into the upper half of the W plane through Schwarz-Christoffel transformation.The approaches of W plane to T plane and T plane to K plane are similar to that from Z plane to W plane.The individual transformations between the planes can be written as follows: where g = ln(R s /R r ).R s and R r are the radius of the stator and rotor, respectively.a and b are the coefficients of w at the corner points of the slot.Finally, the relationship between flux density in the S and K planes is given by: where λ * is the conjugate of the air-gap permeance with λ a and λ b as its real and imaginary parts.
The flux density B k with its real and imaginary parts B r and B θ represents the field solution in the slotless air gap.
There also exist extensive publications about CMM.Lin al. introduced an analytical method for universal motors based on the actual air-gap field distribution from the field solutions by the CMM [119], and presented a generalized analytical field solution for spoke type permanent magnet machines [120].Hanic et al. [95] proposed a novel method for magnetic field analysis based on the CMM and MEC.Alam et al. [121,122] took the slotting effect, winding distribution, armature reaction and saturation effect into account.They presented an improved CMM for magnetic field analysis in surface-mounted permanent magnet motors considering eccentricities.

Virtual Work Principle
The force in the virtual work principle (VMP) is calculated from a spatial derivative of stored energy.The VMP obtains the air gap magnetic field through the electromagnetic theory and the corresponding boundary conditions.Moreover the energy function of the air gap magnetic field is obtained.The partial differentials of energy function in horizontal direction and vertical direction are obtained respectively.By using the method of electromechanical analysis, the coupling dynamics equation is established.These strong nonlinear equations are solved by mathematical transformation.
The magnetic field energy of the air gap space is calculated first and the magnetic tube energy of air gap is: The magnetic field energy of the air gap space is: The electromagnetic force in the x-direction and y-direction can be obtained by the derivation of energy: The VMP has developed for a long time and the existing publications about this method mainly emerged in the late 1990s.Several authors have calculated spatial derivatives of position-dependent air gap reluctances and formulated an analytical expression for force and torque in terms of these derivatives [123,124].Others have implemented VWP discretely by evaluating the total system energy difference and dividing through the spatial difference [125,126].

2D Finite Element Method
The numerical calculation of electromagnetic field is a problem about solution of partial differential equations.For many practical engineering problems, the application of finite element methods (FEM) based on variational principles is often easier than the direct solution of partial differential equations.The Maxwell stress method and virtual work method are usually used for calculating the UMF in the FEM.The FEM divides the continuous field into finite units and then expresses the solution of each unit by the interpolation function which makes it meet the boundary conditions.Then the solution of the continuous field over the entire field is obtained.Numerous scholars studied the FEM for the magnetic field distribution of electrical machines and the FEM is often applied as the contrast or verification for the analytical methods.The FEM can be generally divided into 2D FEM and 3D FEM from the perspective of investigation scope.Moreover, it can be also classified as static FEM and dynamic FEM which is also named time-stepping FEM.
The 2D FEM was adopted to investigate the electromagnetic characteristics of the switched-flux permanent magnet motor [127] and the squirrel-cage induction motor [128].In addition, the time steeping FEM was developed to solve the transient magnetic field in induction machines [129,130].Wang et al. [131] overcame difficulties of hundred stator slots operating under loads and proposed a method for the 2D finite element calculation of UMF in large hydro-generators.Zarko et al. [132] used the FEM to calculate the UMF of a salient-pole synchronous generator in no-load and loaded conditions.Lee et al. [133] adopted the 2D FEM to model magnetic vibration sources in two 100-kW marine fractional-slot interior permanent magnet machines with different pole and slot combinations.

3D Finite Element Method
As illustrated in Figure 12, 2D FEM simplifies the electrical machine system to be a plane problem, while the 3D FEM investigates the issue in the space.Therefore the results obtaining by 3D FEM are more accurate to some extent due to the fact that less assumptions are made.Lee et al. [134] applied the 3D EFM to study the performance of traditional longitudinal flux and transverse flux permanent magnet motors.Chen et al. [135] developed a finite element solution approach for the analysis of the dynamic behavior and balancing effects of an induction motor system.Sibue et al. [136] studied the current density distribution and losses in a large air-gap transformer composed of toe cores and two windings using homogenization and 3D FEM.Ha et al. [137] investigated the coupled mechanical and magnetic forces by the FEM in the transient state solved in a step by step procedure with respect to time.Faiz et al. [129] proposed a time-stepping FEM that identifies mixed eccentricity and overcomes the difficulty of applying FEM to transient behaviors.
The current density is considered as input of 2D FEM which neglects the harmonics of the stator current.The current and rotor angular position can be appropriately considered in the 3D FEM.Generally speaking, the 3D FEM is more time consuming but accurate compared with the 2D FEM.If the size of electrical machines is large or the requirements of calculation accuracy are harsh, the 3D FEM model is needed.
FEM that identifies mixed eccentricity and overcomes the difficulty of applying FEM to transient behaviors.
The current density is considered as input of 2D FEM which neglects the harmonics of the stator current.The current and rotor angular position can be appropriately considered in the 3D FEM.Generally speaking, the 3D FEM is more time consuming but accurate compared with the 2D FEM.If the size of electrical machines is large or the requirements of calculation accuracy are harsh, the 3D FEM model is needed.

Comprehensive Method
The analytical method which can obtain the origin and frequency characteristics of the electromagnetic force is mainly used to explore the mechanism and the solution speed is fast.However the calculation is simplified based on many assumptions and the model is relatively simple.The numerical method can consider more parameters and the accuracy of calculation results is more reliable than the analytical method.However it requires massive time and the computational efficiency for large as well as complex situations is low.Based on these characteristics, some scholars develop a numerical plus analytical calculation method, which takes advantage of the merits of simple analytical method and numerical method.Numerical method is used to solve electromagnetic field distribution for obtaining accurate magnetic flux density.Analytic method is adopted to get the electromagnetic force and the computational efficiency is improved.In addition, experiments are applied to verify the calculation results and this approach is being widely promoted.He et al. [138] combined the FEM with the boundary finite element method to study the electromagnetic, mechanical and acoustic characteristics of permanent magnet DC motors.First, the MMF harmonic of the candidate winding was obtained by analytical method and then each set of harmonics was analyzed by a simplified FEM.Li et al. [139] proposed a semi-analytical method to analyze the eddy current loss of an axial flux permanent magnet motor.Tudorache et al. [140] proposed a hybrid model (numerical and analytical combination) to reduce the cogging torque of a permanent magnet synchronous motor.Compared with the FEM, the computational time was very low and the calculation accuracy remained high.
Based on the advantage combinations of the analytical and numerical methods, Chao [141] proposed a hybrid method to analyze the UMF in hard disk drive spindle motors.Hanic et al. [95] adapted the conformal mapping and magnetic equivalent circuits to calculate the back electromotive force and cogging torque of saturated surface-mounted permanent magnet machines.Sprangers et al. [142] presented a semi-analytical method based on the harmonic modeling technique and analyzed the magneto-static field distribution in the slotted structure of rotating electrical machines.Guo et al. [143] combined the analytical and FEM techniques to predict the air-gap magnetic field distribution of a permanent motor embedded salient pole.

Main Factors Considered in the Calculation
The values of UMF are influenced by many factors.The existing calculation methods are all based on some assumptions, even though the numerical method (which here refers to FEM) can take numerous situations into account.Considering that proposing an exactly perfect model is rather difficult and may be very complex, some key influential factors are needed to paid more attention.This way is efficient in improving the accuracy and will not increase the complexity much at the same time.According to the literature review, the relatively important factors are magnetic saturation, slot and pole combination and load effects.

Magnetic Saturation
In 1918, Rosenberg [144] discovered the effect of magnetic saturation on UMF, which attracted the attention of many researchers.With the increase of excitation current, UMF increases nonlinearly under normal circumstances, but the existence of magnetic saturation limits the infinite increase of electromagnetic force.In addition, the complex geometry structure of rotor and stator tends to cause magnetic field distortion.The magnetic saturation may also be formed in the narrow places wherein slot and pole are located.Generally speaking, the saturation phenomenon is inevitable because of two major factors.One is due to the saturation effects of magnetization characteristic in the ferromagnetic materials.The other is the existence of magnetic flux leakage caused by the distortion in the narrow air-gap space.There exists much evidence supporting the fact of magnetic saturation [145].The magnetic saturation will have a great influence on the UMF and extensive studies have been conducted.Variable degree of saturation [146] and nonlinear magnetic saturation [147] were investigated, respectively.
Calculation methods considering the influence of magnetic saturation effect are emerging in the last decade.Covo [76] used the slope of the magnetization curve to analyze the effect of saturation on the UMF, which was verified by Tenhunen [67] experimentally.Ohishi et al. [148] made further improvements on this method by applying the magnetization curve of the ferromagnetic material and obtained the polynomial relationship between the air gap flux density and the excitation current in the non-air gap portion of the motor.A polynomial equation describing the magnetization curve of the electrical machine was obtained based on the air gap line.Perers et al. [149] studied the effect of magnetic saturation in a hydro-generator eccentric rotor on UMF, which indicated that magnetic saturation significantly influenced the UMF at high voltage and high load.Dorrell [38] proposed a flexible UMF calculation method considering the situation of axial eccentricity and magnetic saturation, and this approach can be applied in the design process effectively.

Slot and Pole Effects
The N-S poles are common in the electrical rotor.The stator slots with the current form winding current.Many calculation methods ignore the influence of poles and slots on the magnetic field distribution.However, the pole and slot change the density of the magnetic field lines and affect the calculation accuracy of UMF.According to the relationship between the number of rotor poles and the stator slots, the combination of poles and slots can be divided into two major categories: one is the integer combination, that is, the ratio of stator slot number to rotor pole number is a positive integer and the other is a fractional combination.More and more motors begin to adopt fractional combination and the systematic research is necessary to investigate its influence on the UMF.On the other hand, the slot skewing may also have a great influence over the harmonics of UMF [150,151].
Some studies showed that harmonics of air gap magnetic flux density generated by the slotted rotating machinery were an important part of UMF [152].The slotting effect will cause an additional magnetic field [153].Furthermore it was found that the vibration modes and frequency components of higher harmonics strongly depend on the combination of the pole-slot number [154,155] and shape of poles [156].Zhu et al. [60] analyzed the influence of different pole/slot combinations on the UMF and found that UMF increases for higher fractional pole-to-slot ratio.In another article, Zhu et al. [157] further proposed the additive effects (p = 3k + 1, p-pole number and k-slot number) and cancelling effects (p = 3k − 1) of different pole/slot combinations on UMF.The slotting effects and air gap eccentricity are often investigated together in recent years.Zarko et al. [117] proposed an analytical method for magnetic field distribution in the slotted air gap of a surface-mounted permanent magnet motor.Bao et al. [158] combined the effects of eccentricity and slotting to conduct magnetic field monitoring in a submersible motor.According to the basic electromagnetic relationship inside the electrical machine, the magnetic field generated by the stator current and the magnetic field generated by the excitation current are superposed to form the composite air-gap magnetic field under the load condition.Many studies calculated the UMF based on the open-circuit current without load [159].This assumption is applicable and reasonable in ordinary cases because the load has little effect on the calculation of UMF.Wang et al. [2] investigated the effect of different load conditions (no load, half rated and rated load) on the UMF and obtained the conclusion that the load had little effect on the electromagnetic force.Although the load current and power factor have changed with load, the UMF has little change compared with no-load [160].The no load cases are widely investigated and the discussions are extensive.For example, Gaussens et al. [43] proposed a general and accurate approach to determine the no-load flux of field-excited flux-switching machines.Hu et al. [161] presented an improved analytical subdomain model for predicting the no-load magnetic field and cogging force.
However, the load in electrical machines exists actually and the effects of load cannot always be ignored.When the accuracy demands for calculations are harsh and the load is high, the load factor should be taken into consideration seriously.Moreover, in the case of saturation, the existence of load may weaken the UMF.Therefore it should be investigated respectively.Perers et al. [149] pointed out that the UMF decreased with increasing load at a given terminal voltage because the inter-pole leakage flux became more pronounced with increasing load.Zhu et al. [60] established a general analytical model in a two-dimensional polar coordinate system and studied the UMF during loading.Studies showed that the additional UMF was large when the electrical load was high.Dorrell et al. [162] pointed out that there was only rotor flux component when the motor was open without load.Flux components were contained in both stator and rotor under load conditions.

The Magnitude and Frequency of the Electromagnetic Force
Electromagnetic force is an electromagnetic attraction that pushes the rotor toward to the stator.As the vibration increases, the air gap length will become smaller and smaller and forms negative feedbacks which further aggravate the vibration.Therefore the UMF is always in the smallest air gap direction.The magnitude of the force is influenced by many factors such as eccentricity, current and winding structure etc. Research have shown that the magnitude of UMF has the nonlinear relationship with these factors and specific situation should be investigated respectively [163].The UMF is rotating if the electrical machine is on operation.Otherwise the force is static.
The common source of UMF is eccentricity.In addition, the static and dynamic eccentricities tend to coexist.The equation describing the frequency components of interest is [22]: where R is the number of rotor slots, p is the number of fundamental pole pairs, f 1 is the fundamental supply frequency, s is the slip, n d is the eccentricity order (n d = 0 in case of static eccentricity and n d = 1 in case of dynamic eccentricity), n is the any positive integer, n sat models magnetic saturation (n sat = 0, 1, 2, . . . ) and k is the order of the stator time harmonics that are present in power supply driving the motor.When neglecting saturation and considering only static eccentricity, the frequency components that are characteristic of a failure agree exactly with the rotor slot harmonics.Furthermore, in the case of saturation, new sidebands around dynamic-eccentricity-characteristic components will appear, according to the combinations of the feasible values of n d and n sat .Classic theory predicts that spatial harmonics of (p ± 1) pole pairs will result from pure dynamic eccentricity.The frequencies of these harmonics will be given by (f 1 ± f r ), where f r is the rotor frequency.Nevertheless, in practical applications, it is likely that these additional fields caused by dynamic eccentricity may induce currents in stator windings because the motor is not completely electrically and magnetically symmetrical.Moreover, the low-frequency components near the fundamental given by: where f r is the rotating frequency.
The magnitude of UMF due to rotor eccentricity is relatively large because of the existence of the direct-current (DC) component in the static eccentric cases and the low fundamental frequency content in dynamic eccentric cases.The magnitude of the force decreases with the increasing of the frequency harmonic.Therefore the magnitudes of the DC component and low frequency content are large and dominate [164].The effects of first several frequency harmonic are often investigated in detail among publications.Arkkio et al. [165] pointed out that the static eccentricity generates an additional force component varying at twice of the supply frequency.Pennacchi et al. [166] investigated the rotating component and twice of it in the UMF for a three-phase generator.Li et al. [167] found that the 3rd harmonic magnetic force plays a major role in the production of the squeaking noise in small permanent magnet DC brush motors.

Measures to Reduce Electromagnetic Vibration
Electromagnetic vibration will cause stability problems and affect the safe and stable operation of the system.Some measures are needed to reduce the vibration of the electromagnetic excitation source.Currently the equalizing winding, damping winding, parallel circuit, magnetic saturation and slot/pole combination are mainly adopted in the electromagnetic vibration suppression.The parallel circuit and damping winding methods are the effective approaches that have been extensively investigated [27].
The parallel circuit is achieved by reconnection of the stator coils groups.Magnetic field harmonics due to rotor eccentricity generate currents circulating in the parallel paths of the rotor and stator windings.These currents equalize the magnetic field distribution in the air gap and hence reduce the resultant UMF.Burakov et al. [100] compared the two approaches (stator parallel paths and rotor windings) reducing UMF in detail.The study found that parallel stator windings can reduce the UMF more effectively than damping windings.Wallin et al. [168] used parallel circuit method to reduce UMF of a synchronous motor.The research shows that the reduction of UMF strongly depends on the relative unbalanced direction of the stator current isolation line.
In the studies of damping winding, Dorrell et al. [169] found that damped windings can significantly reduce the UMF.Wallin et al. [170] studied the effect of damping winding on UMF in an eccentric salient pole synchronous motor.It was found that continuous or discontinuous damping winding produced different damping winding currents, but the influence on UMF was similar.Dorrell [70] analyzed the effect of damped windings in induction machines on reducing the UMF, and thereby reducing bearing wear on the rotor system.
In addition, many other methods are applied to reduce the UMF.For example, Nguyen et al. [171] designed the dual-stator to reduce the UMF for a wound-field flux switching machine.Bi et al. [172] revealed that lead wires can generate severe UMF in permanent magnet synchronous motor, especially at high speed and propose several ways to reduce lead wire asymmetry.Oliveira et al. [173] proposed an equipotential bonding method to change the magnetic attraction force on the circumference of a rotor, and thereby reducing the UMF.

Experimental Study
The publications on theoretical modeling and analysis about UMF are extensive, but experimental data for the measurement of UMF is sparse due to the difficulties in building suitable experimental equipment.The measurement of the air gap field is rather difficulty.Dorrell et al. [145] reviewed the experimental data available about UMF in induction motors and also put forward further data for consideration.
The experiment method of UMF can be divided into direct approach and indirect approach.Considering that the UMF in a real machine is difficult to be measured, the indirect approach is usually adopted [174].For instance, Pennacchi et al. [16] evaluated the UMF effects indirectly by means of the vibrational behavior of the machine.Zarko et al. [132] carried out the measurement and analysis of bearing and shaft vibrations with no-load and loaded conditions.Kim et al. [175] developed an experimental device to measure the axial UMF and verify the simulated results in brushless DC motors.The experimental setup is illustrated in Figure 13   The direct experiment method was studied by Lee et al. [176].As displayed in Figure 14, the rotor is separated from a stator and the stator eccentricity is adjustable.The disassembled rotor is clamped with a load cell and there are no bearings.This experimental device measures the variation of UMF by the load cell on the stationary rotor.The experimental setup is demonstrated as follows: The experimental approaches are usually adopted to conduct some verifications for the analytical or numerical calculation methods.Kumar et al. [85] compared the analytical results with the experimental results based on the back electromotive force (EMF) which is obtained from the airgap field distribution.Kim et al. [128] compared the back EMFs measured in the experiment with simulation results and verified the proposed FEM method.Wang et al. [177] conducted the experiment to validate the theoretical analyses and FE results for the surface mounted permanent- The direct experiment method was studied by Lee et al. [176].As displayed in Figure 14, the rotor is separated from a stator and the stator eccentricity is adjustable.The disassembled rotor is clamped with a load cell and there are no bearings.This experimental device measures the variation of UMF by the load cell on the stationary rotor.The experimental setup is demonstrated as follows:  The direct experiment method was studied by Lee et al. [176].As displayed in Figure 14, the rotor is separated from a stator and the stator eccentricity is adjustable.The disassembled rotor is clamped with a load cell and there are no bearings.This experimental device measures the variation of UMF by the load cell on the stationary rotor.The experimental setup is demonstrated as follows: The experimental approaches are usually adopted to conduct some verifications for the analytical or numerical calculation methods.Kumar et al. [85] compared the analytical results with the experimental results based on the back electromotive force (EMF) which is obtained from the airgap field distribution.Kim et al. [128] compared the back EMFs measured in the experiment with simulation results and verified the proposed FEM method.Wang et al. [177] conducted the experiment to validate the theoretical analyses and FE results for the surface mounted permanent- The experimental approaches are usually adopted to conduct some verifications for the analytical or numerical calculation methods.Kumar et al. [85] compared the analytical results with the experimental results based on the back electromotive force (EMF) which is obtained from the air-gap field distribution.Kim et al. [128] compared the back EMFs measured in the experiment with simulation results and verified the proposed FEM method.Wang et al. [177] conducted the experiment to validate the theoretical analyses and FE results for the surface mounted permanent-magnet machine.

Accurate Calculation of Electromagnetic Excitation
The mechanism of UMF is complicated and further exploration is needed.Moreover the calculation of UMF is still the focus of future research.There is no universally applicable nonlinear analytical expression for computing the UMF so far.Although extensive research about the nonlinear analytical expression of UMF has been conducted, it is based on some assumptions and simplifications.The computational accuracy needs to be improved and the applicability of obtained expressions are limited.The finite element method is more accurate, but tedious steps and time-consuming highlight the limitations in practical applications.Accurate analysis and calculation of UMF under different operating conditions is an important research topic in the design.Therefore, it may be a future research direction to explore the accurate and universal nonlinear model for UMF.

Control of Electromagnetic Vibration
Electromagnetic vibration threatens the safe operation of the entire system.Measures to reduce or eliminate the UMF for the purpose of reducing vibration and noise have been studied.These measures mainly make the electromagnetic field between the stator and stator be as uniform as possible from the aspects of improving the manufacture and installation precision, carrying out multiple balancing checks and arranging the circuit structure rationally.However, as electrical machines move toward high-speed, heavy loads and subtle trends, efficiency and operability will constrain the use of these measures.If a sudden fault caused by the electromagnetic field changes affects the normal operation of the electrical machines, all possible influencing factors need to be checked.This situation cannot meet the long-time and trouble-free design requirements.The frequency components of UMF can be adopted as the basis for fault diagnosis and the fault characteristics can be extracted for targeted dynamic parameters adjustment.At the same time, the dynamic characteristics of the electrical machine can be applied as the monitoring object for fault diagnosis, and thus forming a feedback system.Real-time monitoring and adjustment of the parameters can ensure safe and stable operation of the electrical machine.Therefore, the adoption of electromagnetic vibration characteristics in the vibration control will be a worthwhile study field in the future.

Multiphysics Coupled Modeling
At present, the research of electromagnetic vibration mainly investigates the effect of UMF on solid structures, and ignores the further influence of the structural field on the electromagnetic field.The weak interaction between solid field and electromagnetic field was involved.Future studies will gradually focus on the strong interaction between solid and electromagnetic fields.The end of electrical machines will heat up during operation, which causes variations in the magnetic field.In addition, the mechanical strength and stiffness properties of the device will change.All these factors influence the vibration characteristics.Moreover the air flow of ventilation components, air flow interference during high-speed operation and the external sound field excitation may affect the dynamic characteristics of electrical machines.Coupling vibration between rotating machinery and electromagnetic field is worth exploring.Meanwhile, the rotating structure is affected by the temperature field and surrounding acoustic excitation [178].Due to the flexibility of structure, the multiphysics interaction will produce significant effects on the dynamic properties of structures.

Figure 2 .
Figure 2. Two common cases of parallel eccentricity: (a) static eccentricity and (b) dynamic eccentricity.

Figure 2 .
Figure 2. Two common cases of parallel eccentricity: (a) static eccentricity and (b) dynamic eccentricity.

Figure 2 .
Figure 2. Two common cases of parallel eccentricity: (a) static eccentricity and (b) dynamic eccentricity.

Figure 3 .
Figure 3. Schematic diagram of static and dynamic mixed eccentricity.

Figure 3 .
Figure 3. Schematic diagram of static and dynamic mixed eccentricity.

Figure 5 .
Figure 5. Schematic diagram of the mixed inclined eccentricity.

Figure 5 .
Figure 5. Schematic diagram of the mixed inclined eccentricity.

Figure 5 .
Figure 5. Schematic diagram of the mixed inclined eccentricity.

Figure 6 .
Figure 6.Schematic diagram of the curved eccentricity.Figure 6.Schematic diagram of the curved eccentricity.

Figure 6 .
Figure 6.Schematic diagram of the curved eccentricity.Figure 6.Schematic diagram of the curved eccentricity.

Figure 7 .
Figure 7. Summary of different air-gap eccentricity types.

Figure 7 .
Figure 7. Summary of different air-gap eccentricity types.

Figure 8 .
Figure 8. Calculation method for the unbalanced magnetic force.

Figure 9 .
Figure 9. Schematic diagram for the exact subdomain analysis.

Figure 9 .
Figure 9. Schematic diagram for the exact subdomain analysis.

Figure 10 .
Figure 10.A portion of magnetic equivalent circuit network in electrical machines.Adapted from [97].

Figure 10 .
Figure 10.A portion of magnetic equivalent circuit network in electrical machines.Adapted from [97].

Energies 2018 ,Figure 11 .
Figure 11.Basic steps for obtaining the field solution in a slotted air gap based on CMM.

Figure 11 .
Figure 11.Basic steps for obtaining the field solution in a slotted air gap based on CMM.

Figure 13 .
Figure 13.Experimental setup used to measure the axial UMF of brushless DC motors.Adapted from [64].

Figure 13 .
Figure 13.Experimental setup used to measure the axial UMF of brushless DC motors.Adapted from [64].

Figure 13 .
Figure 13.Experimental setup used to measure the axial UMF of brushless DC motors.Adapted from [64].
as follows: