**Magnetic Field Characteristics and Stator Core Losses of High-Speed Permanent Magnet Synchronous Motors**

#### **Dajun Tao 1,2,\*, Kai Liang Zhou 1, Fei Lv 1, Qingpeng Dou 1, Jianxiao Wu 1,2, Yutian Sun <sup>1</sup> and Jibin Zou <sup>3</sup>**


Received: 4 January 2020; Accepted: 20 January 2020; Published: 22 January 2020

**Abstract:** This study focuses on the core losses in the stator region of high-speed permanent magnet synchronous motors, magnetic field characteristics in the load region, and variations in iron losses caused by changes in these areas. A two-pole 120 kW high-speed permanent magnet synchronous motor is used as the object of study, and a two-dimensional transient electromagnetic field-variable load circuit combined calculation model is established. Based on electromagnetic field theory, the electromagnetic field of the high-speed permanent magnet synchronous motor under multi-load conditions is calculated using the time-stepping finite element method. The magnetic field distribution of the high-speed permanent magnet synchronous motor under a multi-load condition is obtained, and the variations in iron core losses in different parts of the motor under multi-load conditions are further analyzed. The calculation results show that most of the stator iron core losses are dissipated in the stator yoke. The stator yoke iron loss under the no-load condition exceeds 70% of the total stator iron core loss. The stator yoke iron loss under rated operation conditions exceeds 50% of the total stator iron core loss. The stator loss under rated load operation conditions is higher than that under no-load operation. These observations are sufficient to demonstrate that the running status of high-speed motors is closely related to the stator iron losses, which have significance in determining the reasonable yoke structure of high-speed and high-power motors and the cooling methods of motor stators.

**Keywords:** high-speed permanent synchronous motor; magnetic field characteristic; iron loss; stator structure

#### **1. Introduction**

With the increasing market demand for high-efficiency motors for driving high-speed loads, high-speed permanent magnet synchronous motors (HPMSMs) have attracted considerable attention owing to their high power density and high efficiency [1–3]. In recent years, with the development of power electronics technology, inverters have been widely used in motors [4], allowing the motors to operate over a wide range of frequencies. Hu and Gu proposed an adaptive robust three-step control method to eliminate the influence of cogging torque and model uncertainty on the tracking control of a dc motor when its speed varies nonperiodically. This method provides new ideas for various types of motors [5]. However, when the inverter is in operation, many harmonic currents are introduced, and these harmonic currents will affect the internal magnetic field of the motor [6]. The additional losses associated with these harmonics in HPMSMs are greatly increased. This inverter characteristic will reduce the efficiency of a motor, and this reduction is especially pronounced in HPMSMs. Another problem inherent in the addition of harmonic currents is that the eddy current loss generated in the permanent magnets will cause the local temperature of the permanent magnets to rise along with the required increase in power density for high-speed operation that increases the heat load per unit volume. These two sources of additional heat generation present greater difficulties for the heat dissipation of the motor and will even cause local irreversible demagnetization to occur, which will affect the service life and reliability of the motor [7,8].

The losses caused by harmonics when the inverter is operated are called additional losses, and they contain four parts: -1 additional losses in the windings. When the higher order harmonic current passes through the winding, the current density on its cross-section is distributed to the outer surface of the conductor, thus, the equivalent resistance increases and a skin effect is generated. The resulting loss increase is called additional winding loss; -2 additional core iron losses. The additional iron loss is the increase in iron loss caused by the magnetic field generated by the harmonic currents alternating in the stator and rotor cores; -3 eddy current losses of the permanent magnet. The magnetic fields generated by the harmonic currents are not synchronized with the rotor rotation speed, so that the harmonic magnetic field alternates in the permanent magnet, causing eddy current loss to occur; -<sup>4</sup> the eddy current losses of surrounding structures. The magnetic field generated by the harmonic currents will alternate in the surrounding metal structure, causing eddy current loss. High-speed motors are usually driven by high-frequency inverters, and the voltage output by the inverters is rich in higher-order harmonics [9]. The iron loss caused by high-frequency power supplies accounts for a large proportion of the high-speed motor losses. Accurate analysis of the magnetic field characteristics of high-speed motors operating under multiple operating conditions of high-frequency inverter power supplies can provide important support for this study and the analysis of motor iron losses, as well as providing the basis for the accurate calculation of motor core iron losses.

Many scholars, both foreign and domestic, have performed detailed research on motor iron loss [10–20]. Much of this research has been concerned with analytical methods and modulation ratios. Sun Ming's research [10] established an analytical model of the no-load air-gap magnetic field of the axial flux permanent magnet motor, carried out the analytical calculations, and compared the results with finite element calculation results. Some researchers [11,12] analyzed the root cause of rotor eddy current loss, and calculated this loss using the finite element method. The authors in [13] studied the influence of the modulation ratio and carrier ratio on stator losses of permanent magnet synchronous generators by using the two-dimensional field-circuit coupling time-step finite element method. A study [14] deduced the analytical algorithm of the motor core loss under the power supply of a Pulse Width Modulation (PWM) inverter and determined the relationship between the modulation ratio and the iron loss from the perspective of the analytical formula. The results showed that the larger the modulation ratio, the smaller the iron loss. There are also many studies detailing the relationship between fundamental waves and harmonics. The authors in [15] studied the distribution characteristics of the fundamental wave and harmonic iron loss in the stator and rotor core of asynchronous motors, analyzed the waveforms of the magnetic flux density at different positions of the iron core over time during no-load operation, and obtained the different areas of the iron core loss distribution of iron consumption. The authors in [16] performed a harmonic analysis of the magnetomotive force of the motor and used the two-dimensional finite element method to calculate the eddy current losses of the stator and the rotor core and the eddy current losses in the permanent magnet. Yamazaki built a time-step finite element model of a high-speed asynchronous motor and analyzed the iron core loss caused by higher harmonic magnetic fields [17]. Stators and rotors of different materials have specified reference values for the study of iron loss. Many scholars are also studying these various materials and their properties. Denis [18] carried out experiments on a permanent magnet synchronous motor using a nanocrystalline magnetic material as a stator. The results show that compared with an equivalent motor using a conventional non-oriented silicon steel stator core, the nanocrystalline stator

reduced total iron loss to 64% from 75%. Okamoto [19] used a stator core of amorphous magnetic material instead of a stator core of non-oriented steel. The core loss of the motor in this study was reduced by about 50%. The numerical calculation and experimental test data were compared, and the accuracy and reliability of the results was verified. Guo introduced the core loss calculation of the magnetic flux change in a permanent magnet transverse flux machine with a soft magnetic composite stator core and a low carbon steel rotor yoke, based on a modified core loss model and finite element magnetic field analysis. The calculation of the motor core loss is consistent with the measured value on the prototype [20]. After analyzing these literature sources, it can be found that many scholars have performed in-depth research on the calculation models and methods of high-speed permanent magnet motors and high-speed induction motors. However, few studies have detailed calculations of the losses in various areas of the stator. There is still much work to be done to study the magnetic field characteristics of the HPMSMs under heavily loaded conditions and the resulting changes in iron loss.

In this study, a 120 kW high-speed permanent magnet synchronous motor is used as the research object. By establishing a two-dimensional transient electromagnetic field-variable load circuit joint calculation model and using the time-step finite element method, the magnetic field distribution characteristics of the motor under rated operating conditions, no-load conditions, and different transition states from no-load operation to normal conditions can be determined. Additionally, regular research can be performed and iron loss distribution in each structural area of the stator can be quantitatively analyzed. Then, the effects of different load conditions on the losses in each structural area of the stator core can be compared and analyzed, providing a reference for more efficient operation and structural design of HPMSMs in various fields.

#### **2. Establishment of High-Speed Permanent Magnet Synchronous Motor Model**

#### *2.1. Motor Parameters and Physical Model*

The prototype parameters of the 120 kW high-speed permanent magnet synchronous motor studied in this paper are shown in Table 1. This type of motor is mainly used to drive high-speed fans using a direct drive mechanism to achieve efficient operation of the fan system. The physical model of the motor is shown in Figure 1.


**Table 1.** Parameters of prototype motor.

**Figure 1.** Solution domain physical model.

#### *2.2. Solving Equations and Boundary Conditions*

The boundary conditions are also marked in Figure 1. The vector magnetic potential *A* is used to analyze the magnetic field of the motor. *A* has only an *AZ* component. *A* has no *x*- or *y*-axis components and also satisfies the nonlinear Poisson equation. The boundary value problem is as follows:

$$\begin{cases} \frac{\partial}{\partial x} \left( \frac{1}{\mu} \frac{\partial A \underline{\chi}}{\partial x} \right) + \frac{\partial}{\partial y} \left( \frac{1}{\mu} \frac{\partial A \underline{\chi}}{\partial y} \right) = -J\_Z + \sigma \frac{dA \underline{\chi}}{dt} \\\ A Z \big| \overline{\underline{\chi} \underline{\chi}} = 0 \\\ A\_Z \big| \overline{\underline{\chi} \underline{\chi}} = A\_Z \big|\_{\overline{\underline{\chi} \underline{\chi}}} \end{cases} \tag{1}$$

where *AZ* is the *z*-axis component of the vector magnetic potential, *AC*, *AB*, and *CB* are the outer boundary of the stator, *JZ* is the source current density, μ is the permeability of the material, and σ is the conductivity of the material.

#### *2.3. Determination of Stator Core Loss Area*

Figure 2 is a cloud diagram of the magnetic induction intensity distribution when the motor is running at no-load. The magnetic density distribution in each area of the motor structure directly affects the core loss distribution [21]. To thoroughly investigate the core loss distribution in different areas of the motor, the stator core is divided into a stator yoke, a stator tooth root, a stator tooth body, and a stator tooth top according to the small diagram marked by the arrow in Figure 2. Then, the iron consumption in each area can be calculated separately.

**Figure 2.** Magnetic flux density distribution of motor during no-load working condition.

#### *2.4. Determination of Study Location*

In order to systematically study the comprehensive magnetic flux properties of the motor, the following components mush be considered: the magnetic flux density of the 1/3 position of the stator teeth of the motor, the magnetic flux density at the center of the air gap, the flux density in the middle region between the permanent magnet and the rotor edge, the magnetic flux density in the region of the rotor magnetic bridge, the magnetic flux density of the stator yoke, the magnetic flux density of the stator tooth position, and the magnetic flux density of various areas of the stator in multi-load states. Figure 3 shows the seven arcs that correspond to the above areas in order to extract the magnetic density. The stable moment value is extracted in the vector direction of the magnetic density. The positions A, B, C, D, E, and F are centered around the rotor axis in the following arrangement: the length from the position of 1/3 tooth of the stator to the center of the circle, the length from the center of the air gap to the center of the circle, the length from the center of the area between the permanent magnet and the rotor edge to the center of the circle, the length from the center of the rotor magnetic isolation bridge to the center of the circle, and the length from the center of the stator yoke to the center of the circle, respectively. The length from the center of the stator tooth to the center of the circle is the radius where an arc can be drawn with an angle of 180◦ in a turn. Position G is a straight line drawn from the inner

center of the stator to the outer center of the stator. The magnetic density of each of the above locations is derived.

**Figure 3.** Study position of the model.

#### *2.5. Numerical Calculation of Core Loss*

Commonly used iron core loss calculation models include the Steinmetz model and the iron core loss separation model. This study uses the iron loss separation model to analyze the iron loss of high-speed motors.

The iron core loss separation model is also called the constant coefficient trinomial model [15], in which the iron loss is composed of three parts: hysteresis loss, classical eddy current loss, and abnormal eddy current loss. The model is expressed as follows:

$$\begin{split} p\_{ft} &= p\_h + p\_\varepsilon + p\_\varepsilon \\ &= k\_l f B\_m^{-2} + k\_\varepsilon \frac{1}{l} \int \left(\frac{d\mathbf{B}}{dt}\right)^2 dt + k\_\varepsilon \frac{1}{l} \int \left(\frac{d\mathbf{B}}{dt}\right)^{1.5} dt (\mathcal{W}/\text{kg}) \end{split} \tag{2}$$

where *pf e* is core loss, *ph* is hysteresis loss, *pc* is classical eddy current loss, and *pe* is abnormal eddy current loss. *kh* is the hysteresis coefficient, *f* is the alternating frequency of the magnetic flux density, *Bm* is the magnetic density value, *kc* is the classical eddy current coefficient, *ke* is the abnormal eddy current coefficient, and *B* is the magnetic flux density.

In steady state operation, the hysteresis loss mainly depends on two factors, the area surrounded by the hysteresis ring and the alternating frequency of the magnetic flux density. The classical eddy current loss and the abnormal eddy current loss depend on the rate of change of the magnetic flux density. The material properties given by the material manufacturer are obtained using a numerical fitting method. The core loss per unit mass obtained using Formula (2) is multiplied by the core mass to obtain the overall core loss.

#### **3. Magnetic Field Characteristics and Loss Analysis under Multiple Operating Conditions**

#### *3.1. Motor Electromagnetic Field Characteristics at No-Load*

In this paper, a 120 kW HPMSM is used as a test example and the time-step finite element method is used to analyze the magnetic density distribution of the no-load state operation. The magnetic flux density and magnetic field line distribution of the motor at no-load are shown in Figure 4.

**Figure 4.** Motor magnetic field distribution during no-load operation.

The magnetic density waveforms extracted according to the positions A, B, C, and D are shown in Figure 5.

**Figure 5.** *Cont.*

**Figure 5.** Magnetic flux density in each region of the motor during no-load operation.

Analysis can be obtained for the following:


Fourier decomposition is performed on the separated air gap radial magnetic density, and the distribution of the harmonic amplitudes with their orders is shown in Figure 6.

**Figure 6.** Harmonic distribution of air gap tangential flux density.

It can be seen from the analysis of Figure 6 that the harmonic amplitudes of the even harmonics are very small and may be neglected. The odd harmonics from 1 to 13 are separated from the total characteristic waveform. The distribution of the fundamental wave, each odd harmonic wave, and the air gap radial magnetic density in a period is shown in Figure 7.

**Figure 7.** Distribution of each odd harmonic in a period.

The harmonic amplitudes of the 1st to 13th odd waves are extracted and the magnitude of the amplitude of each odd wave is compared to the amplitude of the fundamental wave, as shown in Table 2.

**Table 2.** Magnitude of each harmonic and its proportion in the fundamental.


According to the analysis in Table 2, under the no-load operating conditions, as the harmonic order increases, the proportion of the amplitude of the odd wave to the amplitude of the fundamental wave is generally reduced. Specifically, the amplitude of the third harmonic is the highest compared to the amplitude of the fundamental wave. The remaining odd harmonic amplitudes account for a small proportion of the fundamental amplitude.

#### *3.2. Motor Iron Loss Distribution at No-Load*

Table 3 shows the iron core loss of each area of the stator and rotor under no-load operating conditions.


**Table 3.** Stator core loss distribution of motor during no-load operation.

According to the analysis in Table 3, under no-load operating conditions, the iron core loss of the stator yoke accounts for the largest proportion of the total iron core loss of the motor, 73.43%, and the iron core loss in the stator tooth body and stator tooth root area account for 14.56% and 9.58%, respectively. The stator core top core loss accounts for 2.42% of the total motor core loss and the rotor iron core loss is very small accounting for the smallest proportion of the total motor iron core loss, only 0.03%. Thus, the stator iron core loss accounts for the main part of the total core iron loss of the motor, with a proportion is as high as 99.7%. The distribution of the motor iron core loss under no-load operating conditions is shown in Figure 8.

**Figure 8.** Stator core loss distribution in each region during rated operation.

#### *3.3. Analysis of the Internal Magnetic Field Distribution of a Motor under Multiple Load Conditions*

A 120 kW HPMSM is again used as a test example and the time-step finite element method is used to analyze the magnetic density distribution of the multi-load state operation. Figure 9 shows the different magnetic flux densities in different areas of the motor under multiple load conditions.

**Figure 9.** Magnetic induction intensity distribution of motor in multi-load condition.

Using the analysis position determined in Figure 3, the magnetic density waveforms of position E and position F can be obtained as shown in Figure 10, and a comparative analysis can be found in the following:


**Figure 10.** *Cont.*

**Figure 10.** *Cont.*

**Figure 10.** Magnetic flux waveforms of stator teeth and yoke under four operating conditions.

In order to explain the change of magnetic density more accurately, according to position G, the average value of the magnetic density of the center position of the four areas of the stator core under the following four operating conditions were recorded: no-load, 0.6 times rated load, rated load, and 1.1 times rated load. These results are shown in Table 4.


**Table 4.** Average magnetic density of each area on the stator side under four operating conditions.

According to the analysis in Table 4, the average value of the magnetic flux density in each area of the stator is the largest during no-load operation and the average value of the magnetic flux density in each area of the stator is the smallest when running at 1.1 times the rated load. Under all four operating conditions, the average value of the magnetic density of the tooth body is the largest among the four areas of the stator. When the load is increased, the average value of the magnetic density of each area decreases and the magnetic density of the tooth body area has the largest change compared to the other three positions. Comparing the no-load operating condition with 1.1 times the rated load operating condition, the average value of the stator tooth position magnetic density under no-load operating conditions increases by approximately 0.3 T, and the average magnetic density of the yoke is the smallest. The magnetic density of the yoke area has the smallest change compared to the other

three positions, and the average magnetic flux density of the yoke is relatively unchanged under the four operating conditions.

#### *3.4. Research on the Relationship between Load Condition and Stator Core Loss*

If the waveform when the motor is stable and the iron core loss waveform when stable is extracted, then the average value of the extracted portion can be used as the core loss calculation value. The rotor iron core loss is very small when the motor is running, so we will not study it below. Table 5 shows the iron core loss distribution data for different positions of the stator.

**Table 5.** Stator iron core loss distribution in each region under four kinds of operating conditions.


The relationship between the iron consumption of each area of the motor stator and the load is shown in Figure 11.

**Figure 11.** Stator iron core loss changes in each region under four types of operating conditions.

It can be seen from the analysis shown in Figure 11 that with the increase of the load, the core loss of the stator tooth top, tooth body, tooth root, and yoke parts also increase. Specifically, when the no-load operating condition is increased to 1.1 times the rated load operating condition, the increase of the core losses of the stator tooth root and stator yoke are 336.27% and 125.37%, respectively, which are smaller than those of the stator tooth top and the stator tooth body, which are 539.28% and 494.52%, respectively. From the no-load operating condition to 0.6 times the rated load operating condition, the total increase of the stator core loss is the largest. The increase of the core loss in the four areas of the stator tooth top, tooth body, tooth root, and yoke are 383.67%, 388.61%, 241.6%, and 82.3%, respectively. When the rated load operating condition is changed to 1.1 times the rated load operating condition, the total increase of the stator core loss is the smallest, and the increase of the core loss in the four areas of the stator is the smallest. In particular, the stator tooth top and the tooth root both increase less than 5 W, which is almost unchanged.

Through the research in this paper, the loss distribution ratio of each area of the stator is determined. The main part of the stator iron loss of the motor is the stator yoke iron loss. With the calculation of the motor temperature field, the distribution of the stator heat source can be more accurately understood and will make the temperature calculation of HPMSMs more accurate. Figure 12 shows a flowchart of the motor iron loss-temperature field transition.

**Figure 12.** Flow diagram of motor iron loss-temperature field transition.

#### **4. Conclusions**

In this study, a 120 kW high-speed permanent magnet synchronous motor was used as a test example. A two-dimensional transient electromagnetic field-variable load circuit joint calculation model was established, and the time-step finite element method was used to analyze the magnetic field distribution characteristics and laws of the motor under rated operating conditions, no-load conditions, and transitory operating conditions. The four operating conditions studied were no-load, 0.6 times rated load, rated load, and 1.1 times rated load. Each state was quantitatively analyzed to determine the iron loss distribution in each structural area of the stator, and the impact of different load conditions on the losses in each structural area of the stator core of the motor was then compared and analyzed, leading to the following conclusions:


**Author Contributions:** Conceptualization, K.L.Z. and Y.S.; Methodology, F.L.; Software, Q.D.; Validation, D.T., K.L.Z. and F.L.; Formal Analysis, J.W.; Investigation, J.Z.; Resources, D.T.; Data Curation, K.L.Z.; Writing—Original Draft Preparation, D.T.; Writing—Review & Editing, D.T. and K.L.Z.; Visualization, Y.S.; Supervision, F.L.; Project Administration, D.T. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the National Natural Science Foundation of China (grant nos. 51777048 and 51407050).

**Acknowledgments:** This work was supported in part by the National Natural Science Foundation of China (grant nos. 51777048 and 51407050). The authors would like to thank the anonymous reviewers for their valuable comments and suggestions that strengthened this paper.

**Conflicts of Interest:** All of our authors declare together that there is no conflict of interest.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Review* **Review of Time and Space Harmonics in Multi-Phase Induction Machine**

#### **Vladimir Kindl 1,\*, Radek Cermak 1, Zelmira Ferkova <sup>2</sup> and Bohumil Skala <sup>1</sup>**


Received: 13 December 2019; Accepted: 15 January 2020; Published: 19 January 2020

**Abstract:** Modern multiphase electric machines take advantage of additional degrees of freedom for various purposes, including harmonic current injection to increase torque per ampere. This new approach introduces a non-sinusoidal air gap flux density distribution causing additional technical problems and so the conventional assumptions need to be revised. The paper presents a methodology for synthesis of air gap magnetic field generated by a symmetrically distributed multiphase windings including the rotor field reaction due to the machine's load. The proposed method is suitable either for single-layer or double layer windings and can be adopted either for full-pitched or chorded winding including slots effects. The article analyses the air gap flux density harmonic content and formulates conclusions important to multiphase induction motors. It also discusses effects of time harmonic currents and illustrates the principle of changing number of pole-pairs typical for harmonic currents being injected to increase torque.

**Keywords:** multiphase; induction; motor; space harmonics; time harmonics; injection

#### **1. Introduction**

Over the past few years the modern industry has recorded a huge and intensive development in power electronics and drives and has brought many technical upgrades not only in public transportation. This rapid technological progress results in less industrial energy consumption and improves environmental issues. One of the most frequently discussed topics relating to the traffic environment [1–3] is replacing conventional combustion engine vehicles with fully electric (battery) vehicles (EVs). The initial concept of EV [4,5] comes from the previous experience with hybrid electric vehicles (HEVs), which later evolved into the popular plug-in hybrids (PHEVs). No matter what type of EV is considered, they are mostly powered by standard ac three-phase electric motors [6].

More than 80% of electric vehicles currently use rare earth permanent magnet synchronous motors (PMSM) allowing vehicle manufacturers to increase the efficiency [7,8] compared to traditional induction motors, especially at lower motor speeds. The improved efficiency brings either an increase in transport range or a reduction in battery size (and also the cost) for a given motor weight and vehicle specification. Compared to typical city-driving vehicles, powerful highway vehicles are designed to operate at higher speeds, which somewhat reduces the relative performance increase gained by rare earth magnets and reduces their advantage over induction machines [9]. Moreover, the rising cost and complicated political situation on the market with rare-earth magnets is another very important reason [10] why manufactures start choosing induction motors [11] as main drive units for their electrical vehicles. This tendency can be observed especially when talking about higher performance and luxury vehicles, such as e-buses, locomotives, or vehicles from Audi or Tesla.

The modern e-mobility trends increasingly focus on multiphase variable-speed motor drives [12] since they could provide the transportation with numerous traction and economic advantages [13]. Probably the first mention of a multiphase electric drive dates back to 1969, where the author of [14] proposed the concept of a five-phase inverter-fed induction motor lately extended to a six-phase double-star induction motor, referred in [15,16].

This early interest in multiphase machines was mainly initiated by the possibility of torque ripple reduction, higher reliability and higher fault tolerance. Another strong argument is that for a given motor power, the input power per phase is reduced, resulting in lower demands on the inverter power electronics components. Moreover, the winding I2R losses are inversely proportional to the square of winding distribution factor, hence the higher number of phases may significantly increase the overall motor efficiency.

Maximum theoretical value of winding loss reduction is determined in [17] as 8.8%. Moreover, as shown in [18–22], the multiphase motors may be powered with additional time harmonics injected into the winding to decrease the input current (rms value) while keeping the same torque, which is particularly relevant to the traction battery applications.

The main disadvantage lies in limited slots number available for given stator diameter, the greater number of phases, the lower the number of slots per pole per phase and consequently higher magnitudes of space harmonics. This may be even more problematic in case of outer-rotor machines with narrow-shaped slots.

No matter what number of phases is considered, electric motors are always accompanied by non-linearities and parasitic effects, usually connected with harmonics [23–26], that cause ripple of the input current and the torque, produce noise and increase losses.

These harmonics are often produced by dead times in pulse width modulation (PWM), supply voltage unbalance, magnetic circuit saturation, non-sinusoidal winding distribution, lamination slotting and some other non-linearities and asymmetries.

According to the nature of their origin, we can classify these harmonics as time harmonics and space harmonics [27], while their mutual interaction cannot be neglected [28]. This study describes/recapitulates occurrence and behavior of time/space harmonics in multiphase induction machines.

#### **2. Space Harmonics in the Air Gap Magnetic Field**

Any symmetric *<sup>m</sup>*-phase (*<sup>m</sup>* <sup>∈</sup> <sup>Z</sup>) induction machine has a space displacement between any two successive stator phases equal to 2π/*m*. The stator winding is designed as sinusoidally distributed as possible and is fed with balanced *m*-phase sinusoidal currents. The combined effect is equivalent to having the same winding excited with a constant current and rotating at the stator frequency (rotating field established). Ideally, when the number of stator slots *Q*<sup>1</sup> approaches infinity (*Q*<sup>1</sup> → ∞), and no iron core saturation will appear, the winding forms a sinusoidal magnetic field (or magneto-motive force, mmf) in the air gap δ. However, the practical windings are placed into the finite number of slots and the machines' core always experiences saturation; therefore, the resulting magneto-motive force has rather stepped than sinusoidal curve. This complex curve can be described with Fourier series of mmf waves called space harmonics. The orders of these harmonics are usually marked with symbol ν.

#### *Synthesis of Air Gap Magnetic Field Formed by Symmetrically Distributed Windings*

As shown in [29], a hypothetical single-coil winding, fed by a time-varying sinusoidal current, produces mmf having rectangular waveform according to (1). The equation considers constant air gap permeance independent of angular position (even air gap) and the iron core made of steel having infinite relative permeability.

$$H\_{\mathbf{x}}(a) = \frac{2i}{\pi \delta} \sum\_{\nu=1}^{\infty} \frac{1}{\nu} \sin(\nu \frac{a\_{\mathbf{y}}}{2}) \cos(\nu[a - (\xi + (\mathbf{x} - 1))a\_1]), \ a\_{\mathbf{y}} = \beta \pi \tag{1}$$

In (1), β ∈ < 0; 1 > represents the shortening of the winding coil pitch (chorded winding), *i* is the current content in the slot, α is mechanical angle measured in the air gap periphery, ξ is position of the coil group origin given in slots number and α<sup>1</sup> = 2π/*Q*1.

For symmetric *m*-phase winding distributed in *Q*<sup>1</sup> stator slots, creating 2*p* magnetic poles, we can introduce *q* = *Q*1/2*pm* as the number of slots per phase per pole. Hence, the field produced by respective coils in group of *q* coils can be described using (2). For further analysis, it is reasonable to consider only the basic two-pole winding, therefore we always assume that 2*p* = 2.

$$\begin{array}{l} H\_{1}(\alpha) = \frac{2i}{\pi b} \sum\_{\nu=1}^{\infty} \frac{1}{\nu} \sin(\nu \frac{a\_{\nu}}{2}) \cos(\nu \alpha - \xi \alpha\_{1})\\ H\_{2}(\alpha) = \frac{2i}{\pi b} \sum\_{\nu=1}^{\infty} \frac{1}{\nu} \sin(\nu \frac{a\_{\nu}}{2}) \cos(\nu[\alpha - (\xi + 1)\alpha\_{1}])\\ H\_{3}(\alpha) = \frac{2i}{\pi b} \sum\_{\nu=1}^{\infty} \frac{1}{\nu} \sin(\nu \frac{a\_{\nu}}{2}) \cos(\nu[\alpha - (\xi + 2)\alpha\_{1}])\\ \vdots\\ H\_{q}(\alpha) = \frac{2i}{\pi b} \sum\_{\nu=1}^{\infty} \frac{1}{\nu} \sin(\nu \frac{a\_{\nu}}{2}) \cos(\nu[\alpha - (\xi + (q - 1))\alpha\_{1}]) \end{array} \tag{2}$$

The summation of particular fields (2) gives resulting mmf (3) generated by a group of *q* coils corresponding to the coil group of one stator phase.

$$H\_{\rm grup}(\alpha) = \frac{2i}{\pi \delta} \sum\_{\nu=1}^{\infty} \left[ \frac{1}{\nu} \sin(\nu \frac{\alpha\_y}{2}) \sum\_{k=1}^{q} \left[ \cos(\nu[\alpha - (\xi + (k-1))\alpha\_1]) \right] \right] \tag{3}$$

Simplifying (3) we obtain more useful Equation (4).

$$H\_{\rm group}(\alpha) = \frac{2qi}{\pi \delta} \sum\_{\nu=1}^{\infty} \left[ \frac{1}{\nu} \sin \left( \nu \frac{\alpha\_{\overline{\nu}}}{2} \right) \frac{\sin \left( q \frac{\alpha\_{\overline{1}}}{2} \nu \right)}{q \sin \left( \frac{\alpha\_{\overline{1}}}{2} \nu \right)} \cos \left( \alpha \nu - \nu (q - 1 + 2\xi) \frac{\alpha\_{\overline{1}}}{2} \right) \right] \tag{4}$$

From (4), it is easy to find the resulting mmf waveform generated by any symmetric *m*-phase distributed winding designed with integer *q*. We consider twice the number of mathematical phases *m* = 2*m*, therefore ξ = (*k* − 1)2*q* is substituted for the "plus" phases (A, B, C, D, E, ... ), and ξ = *mq* + (*k* − 1)2*q* is substituted for the "minus" phases (A', B', C', D', E', ... ). The input current with angular frequency ω<sup>1</sup> flowing through the *k*-th stator phase is defined as (5):

$$i\_k(t) = I\_{km} \sin\left(\omega\_1 t - \frac{k-1}{m} \, 2\pi\right) \tag{5}$$

Hence by combination (4) with (5) we obtain (6):

$$H\_{m-\text{phase}}(a) = \frac{4q}{\pi\delta}I\_{\text{km}}\sum\_{\nu=1}^{\infty} \left[ \begin{array}{c} \frac{1}{\nu}\sin\left(\nu\frac{a\_y}{2}\right)\sin\left(\nu\frac{ma\_1\nu}{2}\right)\frac{\sin\left(\frac{a\_1}{2}\nu\right)}{q\sin\left(\frac{a\_1}{2}\nu\right)}\\ \sum\_{k=1}^{m} \left[ \sin\left(\nu\frac{2a+(1-k\cdot 4q)a\_1}{2}\right)\sin\left(a\_1t - \frac{k-1}{m}\cdot 2\pi\right) \right] \end{array} \tag{6}$$

Equation (6) can be further modified into (7) to calculate with given coil turns number per slot *N* carrying the input current *I* √ 2. For single-layer winding, *N*/2 must be used instead of *N*.

$$H\_{\rm{m-phase}}(\alpha) = 4\sqrt{2} \frac{N!q}{\pi \delta} \sum\_{\nu=1}^{\infty} \left[ \begin{array}{c} \frac{1}{\nu} \sin(\nu \frac{a\_{\rm{y}}}{2}) \sin(\nu \frac{\max\_{1}}{2}) \frac{\sin(\frac{a\_{1}}{2}\nu)}{q \sin(\frac{a\_{1}}{2}\nu)}\\ \sum\_{k=1}^{\infty} \left[ \sin(\nu \frac{2\alpha + (1-k\cdot 4q)a\_{1}}{2}) \sin(\omega\_{1}t - \frac{k-1}{m}\cdot 2\pi) \right] \end{array} \right] \tag{7}$$

As an example, for the three-phase winding, (7) results in (8),

$$\begin{cases} H\_{3-ph\infty}(\alpha) = 4\sqrt{2} \frac{Nlq}{\pi\delta} \sum\_{\nu=1}^{\infty} \left[ \frac{1}{\nu} \sin(\nu \frac{a\_{\nu}}{2}) \sin(\nu \frac{m\alpha\_{1}}{2}) \frac{\sin\left(\frac{a\_{1}}{2}\nu\right)}{q \sin\left(\frac{a\_{1}}{2}\nu\right)} \right] \sin\left(\nu \frac{2\alpha + (1-4q)\alpha\_{1}}{2}\right) \sin(\omega\_{1}t) + \\ \quad \quad \quad \quad \sin\left(\nu \frac{2\alpha + (1-8q)\alpha\_{1}}{2}\right) \sin\left(\omega\_{1}t - \frac{2}{3}\pi\right) + \sin\left(\nu \frac{2\alpha + (1-12q)\alpha\_{1}}{2}\right) \sin\left(\omega\_{1}t - \frac{4}{3}\pi\right) ] \end{cases} \tag{8}$$

For the five-phase winding we have (9),

$$\begin{split} H\_{5-ph\infty}(\alpha) &= 4\sqrt{2} \frac{Nlq}{\pi b} \sum\_{\nu=1}^{\infty} \left[ \frac{1}{\nu} \sin(\nu \frac{a\_{\nu}}{2}) \sin(\nu \frac{m a\_{1}}{2}) \frac{\sin(q \frac{a\_{1}}{2} \nu)}{q \sin(\frac{a\_{1}}{2} \nu)} \right] \sin(\nu \frac{2\alpha + (1-4q)a\_{1}}{2}) \sin(\omega\_{1}t) + \\ & \quad \sin(\nu \frac{2\alpha + (1-8q)a\_{1}}{2}) \sin(\omega\_{1}t - \frac{2}{5}\pi) + \sin(\nu \frac{2\alpha + (1-12q)a\_{1}}{2}) \sin(\omega\_{1}t - \frac{4}{5}\pi) + \\ & \quad \sin(\nu \frac{2\alpha + (1-16q)a\_{1}}{2}) \sin(\omega\_{1}t - \frac{6}{5}\pi) + \sin(\nu \frac{2\alpha + (1-20q)a\_{1}}{2}) \sin(\omega\_{1}t - \frac{8}{5}\pi)] \end{split} \tag{9}$$

And the seven-phase winding will generate field according to (10):

$$\begin{split}H\_{7-phase}(\alpha) &= 4\sqrt{2}\frac{Nlq}{\pi\delta}\sum\_{\nu=1}^{\infty}\Bigl[\frac{1}{\nu}\sin(\nu\frac{a\_{\text{F}}}{2})\sin(\nu\frac{ma\_{\text{I}}}{2})\frac{\sin(q\frac{a\_{\text{I}}}{2}\nu)}{q\sin(\frac{a\_{\text{I}}}{2}\nu)}\Big]\sin(\nu\frac{2\alpha+(1-4q)a\_{\text{I}}}{2})\sin(a\_{1}t) + \\ &\sin(\nu\frac{2\alpha+(1-8q)a\_{\text{I}}}{2})\sin(\alpha\_{1}t - \frac{2}{\tau}\pi) + \sin(\nu\frac{2\alpha+(1-12q)a\_{\text{I}}}{2})\sin(a\_{1}t - \frac{4}{\tau}\pi) + \\ &\sin(\nu\frac{2\alpha+(1-16q)a\_{\text{I}}}{2})\sin(\alpha\_{1}t - \frac{6}{\tau}\pi) + \sin(\nu\frac{2\alpha+(1-20q)a\_{\text{I}}}{2})\sin(a\_{1}t - \frac{8}{\tau}\pi) + \\ &\sin(\nu\frac{2\alpha+(1-24q)a\_{\text{I}}}{2})\sin(a\_{1}t - \frac{10}{\tau}\pi) + \sin(\nu\frac{2\alpha+(1-28q)a\_{\text{I}}}{2})\sin(a\_{1}t - \frac{12}{\tau}\pi)]] \end{split} \tag{10}$$

The air gap flux density distribution is then obtained by applying (11).

$$B\_{m-\text{phase}}(\alpha) = \mu\_0 H\_{m-\text{phase}}(\alpha) \tag{11}$$

Previous approach assumes a uniform air gap δ, which makes the analyzed air gap field corresponding to the mmf waveform generated by the winding. The stator and the rotor surfaces are slotted in a practical machine, and therefore, the air gap permeance varies along with the machine periphery and generates additional flux waves. Hence, the air gap appears to be slightly wider than its real mechanical size. The widening of the air gap is traditionally considered via Carter's factor *kc* [30]. As proposed in [29], this can be considered by introducing a fictive air gap (12),

$$\delta(\alpha) = \frac{1}{f\_1(\alpha)} + \frac{1}{f\_2(\alpha)} - \delta\_0 \tag{12}$$

where δ<sup>0</sup> represents the initially assumed (even) air gap, and the functions *f*1(α) and *f*2(α) introduce the stator and the rotor slotting, respectively (13).

$$\begin{aligned} f\_1(\alpha) &= a\_0 - \sum\_{\substack{\nu=1 \\ \nu \text{ odd}}}^{\infty} a\_{\nu} \cos(\nu Q\_1 a), \ a\_0 = \frac{1}{k\_1 \delta\_0} \\ f\_2(\alpha) &= b\_0 - \sum\_{\nu=1}^{\infty} b\_{\nu} \cos(\nu Q\_2 a), \ b\_0 = \frac{1}{k\_2 \delta\_0} \end{aligned} \tag{13}$$

The Carter factors *kc*<sup>1</sup> and *kc*2, important to (13), are determined from given stator/rotor slot pitches *td*<sup>1</sup> and *td*<sup>2</sup> using (14),

$$k\_{c1} = \frac{t\_{d1}}{t\_{d1} - \gamma\_1 \delta\_0}, \ k\_{c2} = \frac{t\_{d2}}{t\_{d2} - \gamma\_2 \delta\_0} \tag{14}$$

where γ1,2 comes from (15).

$$\gamma\_{1,2} = \frac{4}{\pi} \left[ \frac{b\_{01,02}}{2\delta\_0} \operatorname{atan} \left( \frac{b\_{01,02}}{2\delta\_0} \right) - \ln \sqrt{1 + \left( \frac{b\_{01,02}}{2\delta\_0} \right)^2} \right] \tag{15}$$

Parameters *b*<sup>01</sup> and *b*<sup>02</sup> represent the slots opening of the stator and the rotor, respectively. To complete substitution into (13), we only need to calculate the values of *av* and *bv* according to (16),

$$\begin{aligned} a\_{\upsilon} &= \frac{4\beta\_1}{\pi \lambda\_{0} \nu} \Big| \frac{1}{2} + \frac{\left(\frac{b\_{01}}{t\_{d1}} \nu\right)^2}{0.78 - 2\left(\frac{b\_{01}}{t\_{d1}} \nu\right)^2} \Big| \sin\left(1.6\pi \frac{b\_{01}}{t\_{d1}} \nu\right) \\\ b\_{\upsilon} &= \frac{4\beta\_2}{\pi \lambda\_{0} \nu} \Big| \frac{1}{2} + \frac{\left(\frac{b\_{02}}{t\_{d2}} \nu\right)^2}{0.78 - 2\left(\frac{b\_{02}}{t\_{d2}} \nu\right)^2} \Big| \sin\left(1.6\pi \frac{b\_{02}}{t\_{d2}} \nu\right) \end{aligned} \tag{16}$$

where β<sup>1</sup> and β<sup>2</sup> come from (17).

$$\beta\_1 = \frac{1}{2} - \frac{1}{2b\_{01}} \left[ b\_{01} \left( \frac{1}{k\_{c1}} - 1 \right) \right], \\ \beta\_2 = \frac{1}{2} - \frac{1}{2b\_{02}} \left[ b\_{02} \left( \frac{1}{k\_{c2}} - 1 \right) \right] \tag{17}$$

However, by substitution (12) into (7), we find that for some chosen *q* the "tooth-to-tooth" synchronization of both waves is not fully observed, and therefore, it is necessary to introduce a correction factor (angular displacement) for one of the coordinate systems. Hence, Equation (7) should be rewritten into (18).

$$H\_{m-\text{phase}}(a) = 4\sqrt{2} \frac{N!q}{\pi \delta(a)} \sum\_{\nu=1}^{\infty} \left[ \begin{array}{c} \frac{1}{\nu} \sin\left(\nu \frac{a\_0}{2}\right) \sin\left(\nu \frac{m\alpha\_1}{2}\right) \frac{\sin\left(\frac{a\_1}{2}\nu\right)}{q \sin\left(\frac{a\_1}{2}\nu\right)}\\ \sum\_{k=1}^{m} \left[ \sin\left(\nu \frac{2a + (1 - k \cdot 4q)\alpha\_1 - a\_{\text{side}}i}{2} \right) \sin\left(\alpha\_1 t - \frac{k-1}{m} \cdot 2\pi\right) \right] \end{array} \tag{18}$$

For any full-pitch winding with odd *q*, we select α*shi f t* = 2π/*Q*<sup>1</sup> in (17), and analogously, for winding with even *q*, we apply α*shi f t* = 0. In case of chorded winding, we use either α*shi f t* = 2π/*Q*<sup>1</sup> for *q* = 1, 2, 5, 6, 9, 10, ... or α*shi f t* = 0 for *q* = 3, 4, 7, 8, 11, 12, ... . To obtain the flux density distribution, formula *Bm*−*phase*(α) = μ0*Hm*−*phase*(α). must be applied to all equations relating to the air gap magnetic field strength distribution.

For the better understanding, we will analyze the air gap magnetic field for three various "fictive" single-layer (or double-layer) windings, considering slotless (7) and slotted (18) motor geometry. This situation is close to no-load motor operation (zero torque and zero rotor field). First, the flux density generated by a three-phase, full-pitched, winding having *Q*<sup>1</sup> = 30, *Q*<sup>2</sup> = 22, 2*p* = 2 and *q* = 5 is shown in Figure 1. The left-side figure presents the flux density distribution as it depends on the air gap angular position, and the right-side figure shows the resulting frequency spectrum. While the red curve depicts the field considering smooth air gap with no slots present on the stator or the rotor, the blue curve shows the situation for slotted lamination.

Second, the flux density formed by a five-phase, full-pitched, winding having *Q*<sup>1</sup> = 30, *Q*<sup>2</sup> = 22, 2*p* = 2 and *q* = 3 is shown in Figure 2.

Finally, the flux density formed by a seven-phase, full-pitched, winding having *Q*<sup>1</sup> = 28, *Q*<sup>2</sup> = 22, 2*p* = 2 and *q* = 2 is shown in Figure 3.

The comparison between discussed windings shows that the motor having more phases can form a smoother magnetic field (lower content of harmonics) than the less phases motor even if it has similar slots number. Hence, this motor can generate less torque ripple, distinguishes lower THDi and can reduce noise. For air gap flux density spectrum, we can formulate several conclusions generally valid for any distributed multiphase winding having integer *q*.

**Figure 1.** (**a**) Air gap field generated by 3-phase winding; field distribution (**left**), (**b**) frequency spectrum (**right**).

**Figure 2.** (**a**) Air gap field generated by 5-phase winding; field distribution (**left**), (**b**) frequency spectrum (**right**).

**Figure 3.** (**a**) Air gap field generated by 7-phase winding; field distribution (**left**), (**b**) frequency spectrum (**right**).

First, the spectrum includes all harmonics orders calculated from (19). When the operator "+" is used, the relevant harmonic generates mmf travelling the air gap together with the fundamental

harmonic but ν-times slower, and to the contrary, the operator "−" gives harmonics producing waves traveling the air gap in opposite direction (also ν-times slower).

$$\mathbf{v} = 2mc \pm 1, \ c \in \mathbb{Z} \tag{19}$$

Second, the spectrum also includes frequency orders called "step" harmonics obtained from (20). They are present in the spectrum mainly due to the fact that the winding is placed in a finite number of stator slots. The stator "slot" harmonics have therefore the same orders as have the "step" harmonics, and hence, they both overlap in the frequency spectrum. The "step" harmonics have the winding factor of the same size as the fundamental harmonic.

$$\nu\_{1strp} = \nu\_{1slot} = c \frac{Q\_1}{p} \pm 1, \ c \in \mathbb{Z}.\tag{20}$$

Besides the "step" harmonic orders, the air gap flux density includes also the rotor "slot" harmonics (21).

$$\nu\_{2\text{slot}} = c \frac{Q\_2}{p} \pm 1, \ c \in \mathbb{Z} \tag{21}$$

Third, as shown in (18) the winding factor considering straight rotor bars is still given by (22),

$$k\_{\rm uv} = \sin(\nu \pi \frac{\beta}{2}) \frac{\sin\left(\nu \frac{\pi}{m'}\right)}{q \sin\left(\nu \frac{\pi}{m'q}\right)}\tag{22}$$

Hence the magnitude of ν-th harmonics can be calculated from fundamental harmonic using (23).

$$B\_{\nu} = \frac{B\_1}{\nu} \, k\_{w\nu} \tag{23}$$

According to previous results, the harmonic orders of the "step" and the stator "slot" harmonics interact and modify the original order (20) magnitudes. Considering only the interaction between the very first "step" and "slot" harmonic orders calculated from (20) when *c* = 1, then the resulting magnitude of (24),

$$\nu = \frac{Q\_1}{p} - 1\tag{24}$$

Becomes, according to (25),

$$B\_{\nu} = B\_{1stcp} \left[ \frac{a\_1}{2a\_0} \frac{Q\_1 - p}{p} + 1 \right] \tag{25}$$

Moreover, analogously, for (26),

$$\nu = \frac{Q\_1}{p} + 1\tag{26}$$

We obtain (27):

$$B\_{\nu} = B\_{1step} \left[ \frac{a\_1}{2a\_0} \frac{Q\_1 + p}{p} - 1 \right] \tag{27}$$

Into (25) and (27), we substitute from (13) for *a*<sup>0</sup> and from (16) for *a*<sup>1</sup> with applying ν = 1. This gives us a rough estimate of analyzed harmonic amplitudes.

To demonstrate validity of the method, we can analyze (using FEA) the air gap magnetic field of real 3 kW five-phase induction motor corresponding to the case study shown in Figure 2. The motor geometry with magnetic distribution is shown on the left side of Figure 4. The right side of the same figure shows the resulting air gap flux density calculated using FEA and compared to given analytical approach. In order to compare comparable, the FEA considers no-load state (zero torque) and magnetic core composed of linear steel having very high relative permeability (μ*<sup>r</sup>* = 104).

**Figure 4.** (**a**) No-load flux density distribution inside the machine, (**b**) air gap flux density.

For the no-load operation, the proposed analytical method of air gap space harmonics prediction works well, but it may fail when the machine operates under load condition. In this case, the induced voltage generates a current flowing through the rotor, which generates its own magnetic field. This rotor field interacts with the original stator field, which in turn produces torque. As the air gap flux density combines both the stator and the rotor magnetic fields, the resulting curve is deformed as compared to the no-load air gap magnetic field. An example of the situation is shown in Figure 5. The motor (from Figure 4) operates under load (3 kW), and its magnetic field is calculated using nonlinear (steel with μ*<sup>r</sup>* = *f*(*I*)) transient analysis to reach as high accuracy as possible. Figure 5 on the left side represents the instantaneous distribution of the magnetic field in the motor cross-section, and the right frame shows the corresponding air gap flux density curve.

**Figure 5.** (**a**) Flux density distribution inside the machine under load, (**b**) air gap flux density.

The no-load and the loaded operation states are compared graphically in Figure 6 by composing the two curves (taken from Figures 4 and 5) in one graph. The red line (taken from Figure 4) shows the flux density curve corresponding to the no-load operation, and the blue curve (taken from Figure 5) shows the magnetic field corresponding to the operation under load.

As obvious from the right side of Figure 6, the frequency spectrum of the loaded flux density includes harmonics that have not been predicted yet, and which are mainly caused by the rotor field presence and by saturation of the motor magnetic core.

**Figure 6.** (**a**) Comparison between no-load and loaded operational state; field distribution (**b**) from FEA, frequency spectrum (right).

The saturated motor flats the air gap magnetic field according to the BH curve of steel used for the motor construction. As the field is a periodic odd function, the spectrum will contain group of "saturation" harmonics (28), even if the stator and rotor are considered to be slot-less.

$$\nu\_{\text{saturation}} = 2(\mathfrak{c} - 1) + 1, \ \mathfrak{c} \in \mathbb{Z} \tag{28}$$

Although we classify these harmonics as space harmonics, their speed and direction do not follow the previous rules. Analyzing the air gap flux density, we find that the resulting (flattened) mmf travels the air gap with a constant shape (slots are not considered), which means that the harmonics created by the saturation must travel at the same speed and in the same direction as the fundamental wave. As a result, these harmonics contribute to the machine's useful torque.

Since these harmonics travel through the air gap synchronously with the fundamental harmonic, their slip must be the same as slip of the fundamental harmonic. Currents induced into the rotor bars will therefore include harmonics corresponding to (28).

Rotor time harmonics can be calculated either by analytical approach based on methodology described in [31] or by using FEA. The authors of [31] represented a squirrel-cage rotor by a star-connected winding with transformation of the end rings into stars. From pre-calculated saturation harmonics [30,31], they derived the induced voltage per leg of the proposed equivalent diagram and calculated the harmonic currents using transformed bar resistance and reactance.

An example relating to motor in Figure 5 is seen in Figure 7, where the left side represents the time dependency of the current (calculated from FEA) flowing through the bar, and the right side shows its frequency spectrum.

Analogously to (1), we can find the magnetic field (29) generated by the *k*th single rotor bar carrying the current *ik*.

$$B\_{\rm bar-k}(\alpha) = i\_k \frac{\mu\_0}{\pi \delta(\alpha)} \sum\_{\nu=1}^{\infty} \frac{1}{\nu} \sin(\nu[\alpha - (k-1)\alpha\_1])\tag{29}$$

A symmetrically manufactured rotor, having bars evenly distributed along the air gap, generates magnetic field of opposite direction to the field formed by the stator. Considering sinusoidal rotor current distribution, Equation (30) describes the resulting rotor field. Here, *I*2*<sup>m</sup>* represents the bar current magnitude.

$$B\_{\rm cage}(\alpha) = \frac{\mu\_0}{\pi \delta(\alpha)} I\_{2m} \sum\_{k=1}^{Q\_2} \left[ \sin(\omega t - p \frac{k-1}{Q\_2} 2\pi) \sum\_{\nu=1}^{\infty} \frac{1}{\nu} \sin(\nu [\alpha - (k-1)\alpha\_1]) \right] \tag{30}$$

*Energies* **2020**, *13*, 496

As shown in Figure 7, the bar current includes, besides the fundamental harmonic component, also an amount of additional time harmonics deforming the previously assumed sinusoidal current wave. Based on this, (30) should be rewritten in (31) to consider *n* rotor time harmonics.

$$B\_{\rm zagr}(\alpha) = \frac{\mu\_0}{\pi \delta(\alpha)} \sum\_{k=1}^{Q\_2} \left[ \sum\_{\mu=1}^n \left[ I\_{\mu} \sin \left( \mu \left| \alpha t - p \frac{k-1}{Q\_2} 2\pi - q \epsilon\_{i\mu} \right| - \varphi\_{\rm slif} \right) \right] \sum\_{\nu=1}^\infty \frac{1}{\nu} \sin(\nu \left[ \alpha - (k-1)\alpha\_1 \right]) \right] \tag{31}$$

New parameter ϕ*i*<sup>μ</sup> represents the phase shift of μth time harmonic component, and ϕ*shi f t* is the angle measured between the stator and the rotor mmfs obtained analyzing the machines equivalent circuit [32,33]. Final magnetic field curve (see Figure 8) is then given as the difference between the stator and the rotor fields, i.e., *Bm*−*phase*(α) − *Bcage*(α).

**Figure 7.** (**a**) Time dependency of bar current (**left**) and (**b**) its harmonic spectrum (**right**).

**Figure 8.** (**a**) Final air gap flux density for motor operating under load condition; flux density distribution (**left**)—blue line inverted, (**b**) frequency spectrum (**right**).

The left side of Figure 8 compares the flux density waveform obtained using proposed analytical approach (red line) to the one derived from FEA (blue line). To ensure better clarity, the blue line is inverted. The right side of Figure 8 shows comparison between harmonic spectrum corresponding to both waveforms.

The results show good agreement between the proposed method and FEA and is therefore an effective way to air gap flux density prediction in any multiphase induction motor having distributed winding with integer slots number per pole per phase.

#### **3. Time Harmonics in Multiphase Winding**

For common three-phase industrial line-started induction machines, the significant source of non-harmonic voltage is mainly the power supply imbalance. Multiphase motors are usually powered from frequency converters used for easy and energy efficient speed control. The inverters (especially the voltage-types) generate voltage having either rectangular or pulsed shape. Thus, its spectrum contains the amount of harmonics μ dependent on the load and the pulse width modulation settings. Non-harmonic power supply can cause parasitic torque ripples, vibrations, increased noise, increased voltage stress of the insulation system and also can generate higher I2R winding losses due to harmonic currents.

Sometimes, particularly when talking about multiphase motors, the time harmonics are injected into the power supply in order to increase torque per ampere. The extra torque is obtained due to the fact that the flux distribution in the air gap is flattened so that the saturation can be avoided for a wider operational range.

#### *3.1. Harmonics Creating 2p Number of Pole Pairs*

Operational properties of significant time harmonics and their influence on the motor can be analyzed composing a time-varying phasor diagram showing the magnetic field (32).

$$F\_{m-\text{phase }\mu} = F\_{m-\mu} \sum\_{k=1}^{m} \left[ \cos \left( \alpha - \frac{k-1}{m} 2\pi \right) \cos \left( \mu \omega\_1 t - \mu \frac{k-1}{m} 2\pi \right) \right] \tag{32}$$

Sine wave current, represented by the second multiplier in the summation, flowing in each of *m* stationary coils, represented by the first multiplier, produces *m* sine varying magnetic fields perpendicular to the rotation axis. The *m* magnetic fields add as vectors to produce a single rotating magnetic field *Fm*−*phase* <sup>μ</sup>.

For example, the fundamental harmonic component working in three-phase motor generates positive sequence field with amplitude (33).

$$\begin{array}{l} F\_{3-\text{phase }1} = F\_{m1} \frac{1}{2} [\cos(\alpha - \omega\_1 t) + \cos(\alpha + \omega\_1 t)] + F\_{m1} \frac{1}{2} [\cos(\alpha - \omega\_1 t) + \cos(\alpha + \omega\_1 t) + F\_{m1} \frac{1}{2} [\cos(\alpha - \omega\_1 t) + \cos(\alpha + \omega\_1 t - \frac{8}{3}\pi)]]\\ \cos(\alpha + \omega\_1 t - \frac{4}{3}\pi)] + F\_{m1} \frac{1}{2} [\cos(\alpha - \omega\_1 t) + \cos(\alpha + \omega\_1 t - \frac{8}{3}\pi)] = \dots =\\ F\_{m1} \frac{3}{2} \cos(\alpha - \omega\_1 t) \end{array} \tag{33}$$

Similarly, we can derive (34) for the 5th and 7th harmonic components obtaining the negative and the positive sequence waves, respectively.

$$\begin{array}{l} F\_{3-\text{phase }5} = F\_{a5} + F\_{b5} + F\_{c5} = \dots = F\_{m5} \frac{3}{2} \cos(\alpha + 5\omega\_1 t) \\ F\_{3-\text{phase }7} = F\_{d7} + F\_{b7} + F\_{c7} = \dots = F\_{m7} \frac{3}{2} \cos(\alpha - 7\omega\_1 t) .\end{array} \tag{34}$$

Performing analyses for randomly chosen motors and harmonics, e.g., the five-phase motor fed by 5th harmonics and seven-phase motor fed by 7th harmonics, we get (35).

$$\begin{array}{l} F\_{5-\text{phase }5} = F\_{m5} \, \frac{5}{2} \cos(\alpha + 5\omega\_1 t) \\ F\_{7-\text{phase }7} = F\_{m5} \, \frac{7}{2} \cos(\alpha - 7\omega\_1 t) \end{array} \tag{35}$$

Normally, the literature recognizes the harmonic content in the voltage curve of the inverter following rule (36),

$$
\mu = 2mc \pm 1, \ c \in \mathbb{Z} \tag{36}
$$

However, in case of harmonic injection drives, the converter can generate any harmonic we need; therefore, (36) could be rewritten in more general (37).

$$
\mu = m\mathfrak{c} \pm 1, \; \mathfrak{c} \in \mathbb{Z} \tag{37}
$$

Based on the previous analyses, we may conclude that by applying the operator "+" on (37), the resulting harmonics create waves traveling the air gap together with the fundamental one, but their speed is μ-times higher. Moreover, to the contrary, the operator "−" gives harmonics producing waves passing through the air gap again μ-times faster but now in opposite direction. Under perfectly balanced conditions, all these harmonics produce rotating magnetic fields having constant amplitudes (38) and angular speeds in time, so they form circular shaped fields.

$$F\_{m-p\text{horse }\mu} = F\_{m\mu}\frac{m}{2} \tag{38}$$

#### *3.2. Harmonics Creating Higher Number of Pole Pairs than 2p*

For any harmonic order originating from (39), the application of (32) gives always zero value. In addition, all phase currents (related to the analyzed harmonic) are in phase with each other and therefore generate zero sequence (non-rotating) field.

$$
\mu = mc, \ c \in \mathbb{Z} \tag{39}
$$

When injecting subharmonic (μ Z), the resulting magnetic field changes the amplitude and angular speed in time, which deforms originally circular field into elliptical. In extreme cases (40), the field takes the pulsating form having zero translation speed against the air gap and generates losses, produces noise and causes homopolar magnetic saturation.

$$
\mu = \frac{m}{2} + m(c-1), \ c \in \mathbb{Z}.\tag{40}
$$

The situation becomes more complicated when studying integer harmonic orders that are not predicted, neither by (37) nor (39). Substituting into (32) we get zero values even though they are harmonics generating fields with different phase shift relative to each other. Therefore, they produce magnetic field passing through the air gap with nonzero speed. For five-phase motor, we find those harmonics from (41),

$$
\mu = m(c - 1) \neq \mathbf{2}, \ c \in \mathbb{Z} \tag{41}
$$

However, the seven-phase motor includes, besides (41), also (42).

$$
\mu = m(c - 1) \pm 3, \ c \in \mathbb{Z} \tag{42}
$$

Hence, in case of nine-phase motor we must apply (41) to (43).

$$
\mu = m(c - 1) \neq 4, \ c \in \mathbb{Z} \tag{43}
$$

All these harmonics produce a magnetic field with higher number of pole pairs than generated by the fundamental harmonic. This is due to the specific "time" re-assembling of the stator phases (μ-multiplication of individual phase shifts) while keeping the same winding mechanical distribution. This in turn forms a new winding with completely new Görges diagram [30] producing higher number of pole-pairs. Figure 9 shows the example of the 3rd harmonic (magnitude equals to the fundamental) injected into the five-phase motor discussed previously in Figure 5 (*Q*<sup>1</sup> = 30, *Q*<sup>2</sup> = 22, 2*p* = 2, *q* = 3).

As it forms 3 × 2*p* number of poles, the speed must be the same as that developed by the fundamental harmonic. Similar behavior can be observed also for the 3rd and 5th current harmonic in seven-phase motor and for the 3rd, 5th and 7th harmonic in nine-phase motor (see left side of Figure 10).

**Figure 9.** (**a**) Stator phases re-assembling, (**b**) mmf for the fundamental and the 3rd harmonics.

**Figure 10.** Time harmonics acting in multiphase winding: (**a**) harmonics suitable for injection, (**b**) harmonics inappropriate for injecting.

Amplitudes of individual harmonics could be derived from new Görges diagram. As these amplitudes are defined by half ampere-turns corresponding to one magnetic pole, we may calculate them from the ampere-turns that belong to one pole pitch. These specific "pole-creating" harmonics generate in the air gap space harmonics having the same order (μ = ν), so we can consider the pole pitch to be μ-times shorter as compared to the fundamental. Hence, only (44) coil turns contribute to the final value of individual harmonic amplitude.

$$q\frac{m}{2\mu}N\tag{44}$$

With multiplication (44) by the mean value of the stator current *I*μ we obtain the resulting mmf from (45).

$$H\_{\mu-\text{mag}} = \frac{q}{\delta\_0} \frac{m}{\mu} \mathcal{N} \frac{I\_\mu}{\pi} \sqrt{2}, \ \mu = \nu \tag{45}$$

Considering the winding factor of μ-th harmonic component (still we assume μ = ν), the magnitude of the air gap flux density is then (46).

$$B\_{\mu-\text{mag}} = \mu\_0 \frac{q}{\delta\_0} \frac{m}{\pi} N \frac{k\_{\text{uv}}}{\mu} I\_{\mu} \sqrt{2}, \ \mu = \nu \tag{46}$$

As demonstrated in the right side of the figure, other harmonics (excluding (37) and (39)) may create either circular or the elliptical fields traveling in different directions and speeds.

Table 1 shows the same overview but extended by nine-phase machines. Column "poles" gives number of magnetic poles created by the particular harmonic component, column "sequence" shows the traveling direction (or field character) and column "speed" specifies the velocity relative to the fundamental wave. The field character shows whether the mmf travels against the air gap ("+" or "−") or pulsates ("puls.") with zero speed.


**Table 1.** Time harmonics effect overview.

Following this, we can construct Figures 11–13 to give graphical demonstration of the harmonics' operational influence.

**Figure 11.** Time harmonics effect demonstration in three-phase machines.

**Figure 13.** Time harmonics effect demonstration in seven-phase machines.

As obvious, there is an amount of harmonics that can create magnetic field traveling the air gap with the same speed and direction that travels the fundamental wave (color marked). These harmonics are usually being injected in order to increase the torque.

#### **4. Conclusions**

From the research results, it is evident that a motor having a higher number of phases can produce a more sinusoidal magnetic field than the motor having a smaller number of phases, even if it has a similar number of slots. Hence, by increasing the number of phases, the torque harmonics decrease. For distributed multiphase winding with integer *q*, we can conclude following rules.

First, the spectrum of the air gap flux density includes harmonic orders obtained from (19). These harmonics generate fields traveling the air gap in a different direction (according to operator used) and speed (ν-times lower than fundamental).

Second, as the windings are physically placed in slots, the spectrum includes significant orders called "step" harmonics (20) having the same winding factor as calculated for fundamental harmonic. The real "slot" harmonics are obtained also from (20), and therefore, their orders overlap (interact) with "step" harmonics. As a consequence, this interaction modifies original magnitudes of (20), and resulting air gap flux density is influenced.

Third, as shown in (18), the winding factor considering straight rotor bars is still given by (22), hence the magnitude of ν-th harmonics can be calculated from fundamental harmonic using (23).

Fourth, the injection of time harmonics can produce magnetic field having even more than 2*p* magnetic pole-pairs. For any harmonic order originating from (39), the application of (32) produces zero value, therefore zero sequence (non-rotating) field.

Fifth, considering subharmonics (μ Z), the resulting magnetic field deforms its originally circular shape into an elliptical one. In extreme case (40), the field pulsates in only one axis forming a kind of homopolar field. When studying integer harmonic orders that are not predicted, neither by (37) nor (39), the substitution into (32) gives zero values even though they are harmonics establishing field waves passing through the air gap win non-zero speed. All these harmonics produce magnetic field having more magnetic pole-pairs than 2*p*. Some of them (e.g., 3rd, 5th and 7th harmonic in nine-phase motor) can travel the air gap with the same speed as developed by the fundamental harmonic and may be therefore injected in order to increase the torque.

**Author Contributions:** Conceptualization, V.K. and R.C.; methodology, V.K.; software, V.K.; finite element analyses, R.C. and Z.F.; formal analysis, V.K.; writing—original draft preparation, V.K.; writing—review and editing, V.K., R.C., Z.F. and B.S.; visualization, V.K.; supervision, Z.F. and B.S. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Ministry of Education, Youth and Sports of the Czech Republic under the project OP VVV Electrical Engineering Technologies with High-Level of Embedded Intelligence CZ.02.1.01/0.0/0.0/18\_069/0009855 and by funding program of the University ofWest Bohemia number SGS-2018-009. This research was also funded by Slovak Research and Development Agency under the contract No. APVV-18-0436.

**Conflicts of Interest:** The authors declare no conflicts of interest.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Acoustic Noise Computation of Electrical Motors Using the Boundary Element Method**

#### **Sabin Sathyan 1,\*, Ugur Aydin <sup>1</sup> and Anouar Belahcen 1,2**


Received: 11 October 2019; Accepted: 31 December 2019; Published: 3 January 2020

**Abstract:** This paper presents a numerical method and computational results for acoustic noise of electromagnetic origin generated by an induction motor. The computation of noise incorporates three levels of numerical calculation steps, combining both the finite element method and boundary element method. The role of magnetic forces in the production of acoustic noise is established in the paper by showing the magneto-mechanical and vibro-acoustic pathway of energy. The conversion of electrical energy into acoustic energy in an electrical motor through electromagnetic, mechanical, or acoustic platforms is illustrated through numerical computations of magnetic forces, mechanical deformation, and acoustic noise. The magnetic forces were computed through 2D electromagnetic finite element simulation, and the deformation of the stator due to these forces was calculated using 3D structural finite element simulation. Finally, boundary element-based computation was employed to calculate the sound pressure and sound power level in decibels. The use of the boundary element method instead of the finite element method in acoustic computation reduces the computational cost because, unlike finite element analysis, the boundary element approach does not require heavy meshing to model the air surrounding the motor.

**Keywords:** acoustics; boundary element method; electric machines; finite element method; induction motors; magneto-mechanics; modeling; noise; vibro-acoustics

#### **1. Introduction**

The acoustic noise in electric motors is a phenomenon of a complex nature and origin. The first kind, electromagnetic vibration and noise, is produced by magnetic forces, magnetostrictive expansion of the core laminations, eccentricity, phase unbalance, slot openings, and magnetic saturation. The second cause of noise is mechanical and is associated with mechanical assembly, in particular the bearings. The third major group is aerodynamic noise, which is due to the flow of ventilating air through or over the motor. These three sources are illustrated in Figure 1. A detailed review on the different forms of vibration and noise in electrical motors can be found in a paper by Vijayraghavan et al. [1]. One form of energy conversion happening in an electrical motor is from electrical energy to acoustic energy. The supply current interacts with the material to produce a magnetic field, which in turn produces magnetic forces. These forces excite the stator core and frame in the corresponding frequency range and produce mechanical vibrations. As a consequence of vibrations, the surface of the stator yoke and frame deforms with frequencies corresponding to the frequencies of forces. These stator and frame vibrations cause the surrounding medium of air to excite and vibrate and finally generate acoustic pressure variations (and thereby noise).

**Figure 1.** Generation of noise of different origins in a rotating electrical machine.

There have been various studies in the literature pertaining to the analysis of the vibrations and noise of electrical motors. Early stages of noise studies used analytical models and combined numerical and analytical methods. Belmans et al. carried out studies on the analytical formulation of acoustic noise in induction motors and authenticated the findings with experimental results [2]. They employed the rotating field theory with the Maxwell tensor method for calculating the frequency components produced by a motor that connect it to the airgap flux density time harmonics produced by the supply. The inference from their study was that high noise levels may be anticipated when one of the frequencies of the electromagnetically excited forces becomes the same as a natural frequency of the stator. They also developed a computerized model using finite element calculations and modal analysis that predicts the frequency components expected in the audible noise of a three-phase induction motor [3]. Cameron et al. have done measurement-based studies on vibrations and noise on reluctance motors and established that the stator deformation due to radial magnetic forces is the leading electromagnetic cause of noise [4]. Besnerais et al. demonstrated the impact of a Pulse Width Modulation (PWM) supply and switching frequencies on the magnetic noise of induction machines using analytical models [5]. Their approach was based on mechanical and acoustic 2D-ring stator models to compute the effect of winding space harmonics and PWM time harmonics in noise production [6]. Later on, Besnerais developed a multiphysical simulation tool for fast calculation of acoustic noise based on analytical and semi-analytical methods [7]. Their platform incorporated different models such as a permeance/magneto motive force (mmf) model, a subdomain model, and a finite element model. The semi-analytical models proved to be faster than the fully finite element models and had the same accuracy level, according to their evaluations. Devillers et al. studied the effect of tangential magnetic forces on vibrations and acoustic noise using a fast subdomain method to calculate Maxwell stress distribution and an electromagnetic vibration synthesis technique [8]. The same team later developed an experimental benchmark set-up for magnetic noise and a vibration analysis of electrical machines [9]. Fakam et al. coupled finite element structural analysis with an analytical tool to compute and compare the electromagnetic noise between surface permanent magnet and interior permanent magnet rotor topologies of a synchronous machine [10]. The same approach of a combined structural finite element method (FEM) and analytical methods was formulated by Islam et al. for computing sound power levels in synchronous motors [11].

With the improvement in numerical computational tools, the use of the boundary element method (BEM) or a combined FEM–BEM in noise computations became widespread. Moreover, these methods can give more accurate results in acoustic calculations. Juhl et al. developed a numerical method based on the BEM to calculate acoustic noise using the Helmholtz integral equation [12]. They created BEM-based numerical software for calculating sound fields on the exterior of bodies of three-dimensional shape or axisymmetric geometries [13]. Wang et al. used the BEM for computing sound power radiating from induction motors and a coupled structural FEM and acoustic BEM in their work [14]. Herrin et al. formulated a high-frequency BEM and compared it to the Rayleigh approximations method. They concluded that the high-frequency BEM is the more robust method [15]. Roivanen has employed different methods such as the BEM, high-frequency BEM, and plate approximation method combined with broad analytical, numerical, and experimental studies for calculating and comparing the sound power levels of electrical motors [16]. Neves et al. [17] presented a study on the coupling between magnetic forces, vibrations, and noise in a switched reluctance motor using FEM and BEM and gathered relatively good results. However, there were some mismatches between measured and simulated results that have been attributed to noise of aerodynamic origin. Furlan et al. [18] followed the same approach of a joint FEM–BEM-based analysis of a permanent magnet Direct Current (DC) electric motor. They simulated all three models, including an electromagnetic analysis in 3D, which could increase the memory requirements and computation time. In this paper, an electromagnetic simulation was done in 2D, and magnetic forces were taken from a 2D model and put into a 3D structural model. Schlensok et al. [19] presented an acoustic simulation of an induction machine with a squirrel-cage rotor, where 3D models were employed for electromagnetic, structural, and acoustic simulations. The electromagnetically excited structure- and air-borne noise of the motor was described in detail in their study. Deng et al. [20] investigated noise in an axial flux permanent magnet motor using electromagnetic and structural FEM and acoustic BEM; however, their approach also included a time-consuming complete 3D scheme. Järvenpää et al. [21] proposed a fast boundary element simulation of noise, where they imposed the surface velocity of a structural FEM model as a source of a fast BEM. This method facilitates the efficient modeling of large acoustic problems, although the computational cost is relatively high in this approach. Besides the computation of sound pressure in pascals and sound levels in dB surrounding the motor (shown in some previous studies), this paper also presents a far-field sound level calculation, which portrays the directivity of sound around the motor.

Inspired by the studies and findings in the literature, this paper investigated the aspect of vibro-acoustics in an electrical motor and formulated a practical and effective approach for acoustic noise calculation. The computational methodology and results successfully present a numerical technique for computing acoustic noise generated by an induction motor. The FEM- and BEM-based model describes how various quantities in different domains in an electrical motor can be calculated and used for acoustic analysis. An extensive numerical analysis of the magneto-mechanical and vibro-acoustic characteristics of a high-speed induction motor used for industrial applications is an original element of this study. Moreover, this paper paves a vivid path for researchers by detailing and clarifying the intricacies related to the preparation of a vibro-acoustic model of an electric motor by explaining both the theoretical and numerical implementation facets of the model. Although many studies have been carried out in this field of research, this paper brings an original contribution through a fully numerical analysis, the implementation of which is explained in detail.

#### **2. Computational Methodology**

This section explains our computational methods, including the main equations and numerical simulation stages employed for calculating the magnetic, mechanical, and acoustic quantities of the motor. Numerical modeling of noise generation in an electrical motor involves three models: first, modeling of the electromagnetic forces; then, modeling of the structural deformation and vibration behavior; and finally, modeling of the consequent acoustic response of the motor, as depicted in Figure 2.

**Figure 2.** Numerical simulation stages.

#### *2.1. Electromagnetics*

The Maxwell equations of the magnetic field problem are solved numerically using 2D FEM. The Maxwell stress tensor method gives the magnetic torque exerted on a ferromagnetic region by integrating the magnetic stress over a surface around it [22]. Equation (1) is

$$T\_{\rm env} = \oint\_{S} r \times \mathbf{r} \cdot d\mathbf{S} = \oint\_{S} r \times \left\{ \frac{1}{\mu\_{0}} (\mathbf{B} \cdot \mathbf{n}) \mathbf{B} - \frac{1}{2\mu\_{0}} \mathbf{B}^{2} \mathbf{n} \right\} dS\_{\prime} \tag{1}$$

where τ is the Maxwell stress tensor, *r* is a vector representing the distance from the integration point to the torque axis, *B* is the magnetic flux density, μ<sup>0</sup> is the permeability of the vacuum, and *n* is the normal unit vector of the integration surface *dS*. In an electrical machine, the integration surface is chosen as the outer surface of the rotor or any cylindrical surface in the airgap.

The magnetic force can be calculated from the Maxwell stress tensor by the volume integral . *<sup>V</sup>* ∇.τ*dV*. This volume integral can be reduced to the closed surface integral over a surface *S*, and the force formula becomes [23]

$$F = \oint\_{S} \left\{ \frac{1}{\mu\_0} (\mathbf{B} \cdot \mathbf{n}) \mathbf{B} - \frac{1}{2\mu\_0} \mathbf{B}^2 \mathbf{n} \right\} dS = \oint\_{S} \left( \frac{1}{2\mu\_0} (\mathbf{B}\_n^{-2} - \mathbf{B}\_t^{-2}) \mathbf{n} + \frac{1}{\mu\_0} \mathbf{B}\_n \mathbf{B}\_t \mathbf{t} \right) dS,\tag{2}$$

where *t* represents the outward unit vector tangential to the differential surface *dS*. The quantity inside the integral is usually interpreted as surface force density or traction. In this study, we follow the same interpretation and compute this force density on a surface located in the inner radius of the stator.

#### *2.2. Structural Mechanics*

Magnetic forces are the excitation parameters for the structural simulation of the stator core of the motor. The forces are fed to the stator elastic model as a body load in the simulation. The elastic model is represented by

$$
\rho \frac{\partial^2 d}{\partial t^2} - \nabla \cdot \boldsymbol{\tau} = \boldsymbol{f}\_{\prime} \tag{3}
$$

where ρ is the mass density, *d* is the vector of displacements, and *f* is the given volume force. After computing the displacements, a discrete Fourier transformation is performed using fast-Fourier transform (FFT), where the time-dependent solution is transformed from times to frequencies in the frequency domain.

#### *2.3. Acoustics*

The BEM used in this study is based on the direct collocation method [13], which deals directly with acoustic variables (sound pressure and particle velocity) and boundary conditions. The multiphysics coupling that combines FEM-and BEM-based physics is employed for coupling the results of solid mechanics finite element physics to acoustics boundary element physics. In the case of an electric motor, the vibrating stator boundary can be used as the acoustics FEM–BEM boundary to couple the acceleration from FEM computation to the BEM interface. This approach allows for modeling in an FEM–BEM framework using the strength of each formulation effectively. Acoustics physics solves the Helmholtz equation for constant-valued material properties and uses the pressure as the dependent variable. The wave equation can be solved in the frequency domain for one frequency at a time. The governing Helmholtz equation for a boundary element interface is given by

$$-\frac{1}{\rho\_c}\nabla^2 p\_t - \frac{k\_{c\eta}^{\quad2}}{\rho\_c} p\_t = 0,\tag{4}$$

$$\left(k\_{c\eta}\right)^2 = \left(\frac{\alpha}{c\_c}\right)^2 \text{ and } p\_l = p + p\_b\tag{5}$$

where *pt* is the total acoustic pressure, *Kcq* is the wave number, ρ*<sup>c</sup>* is the density, ω is the angular frequency, and *cc* is the speed of sound in air.

The governing equations and boundary conditions are formulated using the total pressure *pt* with a scattered-field formulation. In the presence of a background pressure field defining a background pressure wave *pb*, the total acoustic pressure *pt* is the sum of the pressure solved for *p*, which is then equal to the scattered pressure *ps* and the background pressure wave. The equations then contain information about both the scattered field and the background pressure field.

The benefit of the boundary element method is that only boundaries need to be meshed, and the degrees of freedom (DOFs) solved for are restricted to the boundaries. This is beneficial for handling complex geometries, as it introduces easiness in numerical computations. However, the BEM procedure results in fully populated or dense matrices, compared to the sparse system matrices in the FEM. Hence, the BEM is more expensive in terms of memory requirements per DOF than the FEM is, but it has fewer DOFs. Assembling and solving these can be very demanding. In this context, when solving acoustic models of small and medium size, the FEM will often be faster than solving the same problem with the BEM. This could be interpreted as one limitation of the BEM in small-sized computational models. However, in the case of the FEM, when the geometries are complex or two structures are far apart, large air domains need to be meshed. This is costly in terms of the numerical computational facet, as the frequency is increased.

The acoustic pressure computation problem involves solving for small acoustic pressure variations *p* in the surrounding medium of a sound source on top of the stationary background pressure *p*0. In mathematical terms, this can be interpreted as a linearization of the dependent variables around the stationary quiescent values. The fluid flow problems in a compressible lossless fluid can be analyzed using the three governing equations, viz., the mass conservation equation or the continuity equation, the momentum conservation equation or Euler's equation, and the energy equation or the entropy equation. They are given by

$$\frac{\partial \rho}{\partial t} + \nabla \cdot (\rho u) = \mathbf{M}\_{\prime} \tag{6}$$

$$\frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u} \cdot \nabla)\mathbf{u} = -\frac{1}{\rho} \nabla p + \mathbf{F},\tag{7}$$

$$\frac{\partial \mathbf{s}}{\partial t} + \nabla \cdot (s\mathbf{u}) = \mathbf{0},\tag{8}$$

where ρ is the total mass density, *p* is the total pressure, *u* is the velocity field, *s* is the entropy, and *M* and *F* are the possible source terms representing the body forces, if any. In conventional pressure acoustics scenarios, all thermodynamic processes are assumed to be isentropic in nature, which means the processes are both reversible and adiabatic. The small-parameter expansion is executed on a stationary fluid (*u*<sup>0</sup> <sup>=</sup> 0) of density <sup>ρ</sup><sup>0</sup> (kg/m3) and at pressure *<sup>p</sup>*<sup>0</sup> (Pa) such that *<sup>p</sup>* <sup>=</sup> *<sup>p</sup>*<sup>0</sup> <sup>+</sup> *<sup>p</sup>*<sup>1</sup> with *<sup>p</sup>*<sup>1</sup> *<sup>p</sup>*0, <sup>ρ</sup> <sup>=</sup> <sup>ρ</sup><sup>0</sup> <sup>+</sup> <sup>ρ</sup><sup>1</sup> with <sup>ρ</sup><sup>1</sup> <sup>ρ</sup>0, *<sup>u</sup>* <sup>=</sup> <sup>0</sup> <sup>+</sup> *<sup>u</sup>*<sup>1</sup> with *<sup>u</sup>*<sup>1</sup> *<sup>c</sup>*, and *<sup>s</sup>* <sup>=</sup> *<sup>s</sup>*<sup>0</sup> <sup>+</sup> *<sup>s</sup>*1.The small acoustic variations are represented by the variables with subscript 1. Inserting these values in the governing equations gives

$$\frac{\partial \rho\_1}{\partial t} + \nabla \cdot (\rho\_0 u\_1) = \mathbf{M},\tag{9}$$

$$\frac{\partial u\_1}{\partial t} = -\frac{1}{\rho\_0} \nabla p\_1 + F\_\prime \tag{10}$$

$$\frac{\partial p\_1}{\partial t} = \mathfrak{c}\_s^2 (\frac{\partial \rho\_1}{\partial t} + \mathfrak{u}\_1 \cdot \nabla \rho\_0)\_\prime \tag{11}$$

where *cs* is the isentropic speed of sound. The pressure time differential in the last equation is derived from the entropy equation. If the material parameters are constant, the last equation reduces to

$$
\rho\_1 = c\_s^2 \rho\_1. \tag{12}
$$

This expression of acoustic pressure gives a condition that needs to be fulfilled for the linear acoustic equations to hold:

$$\left|p\_1\right| \ll c\_s^2 \rho\_0. \tag{13}$$

Finally, the transient wave equations for pressure waves in a lossless medium can be obtained by rearranging Equations (9)–(11) and dropping the subscripts:

$$\frac{1}{\rho c^2} \frac{\partial^2 p}{\partial t^2} + \nabla \cdot \left[ -\frac{1}{\rho} (\nabla p - q\_d) \right] = Q\_{\text{m}} \tag{14}$$

where the source term *Qm* is a monopole domain source corresponding to a mass source and *qd* is a dipole domain source representing a domain force source. The speed of sound (*c*) and the density (ρ) may in general be space-dependent. The combination term ρ*c*<sup>2</sup> is called the adiabatic bulk modulus (*Ks*) with unit Pa, which is related to the adiabatic compressibility coefficient β*<sup>s</sup>* = 1/*Ks*.

In the frequency domain, the Helmholtz equation can be written as

$$\nabla \cdot \left[ -\frac{1}{\rho\_c} (\nabla p\_t - q\_d) \right] - \frac{k\_{c\eta} \, ^2 p\_l}{\rho\_c} = Q\_{\text{fl}}.\tag{15}$$

Acoustics problems mostly encompass simple harmonic waves such as sinusoidal waves. In numerical computations, to model acoustic–structure interactions, a structural analysis can be coupled to acoustics by imposing the acceleration as a source in the boundaries of the structure in the form of normal acceleration, specified as

$$-\mathbf{n} \cdot \left[ -\frac{1}{\rho\_c} (\nabla p\_t - q\_d) \right] = -\mathbf{n} \cdot a\_{0\prime} \tag{16}$$

where *a*<sup>0</sup> is the normal acceleration and *qd* is the external force term.

#### *2.4. Acoustic Pressure and Audible Sound*

Sound is measured by changes in air pressure. The louder a sound is, the larger the change in air pressure is. The change here is the change from normal atmospheric pressure or reference pressure to the pressure disturbance caused by the sound. Sound pressure is measured in the unit "pascals". A pascal (Pa) is equal to the force of one newton per square meter. The smallest sound pressure a human ear can hear is 20 μPa, which corresponds to zero dB. The sound pressure level (SPL) in dB can be calculated by

$$SPL = 20\log\_{10}\left(\frac{p}{p\_{ref}}\right) \text{dB},\tag{17}$$

where *pre f* is the reference pressure 20 μPa in the case of audible sound calculations.

#### **3. Results**

The results of the numerical simulations are presented in this section. The simulations were done using COMSOL multiphysics software [24]. The specifications of the solid rotor induction motor are given in Table 1. The noise computation of an electric motor starts from the electromagnetic field computations, where the Maxwell equations are solved using finite element analysis. From the electromagnetic analysis, the forces of electromagnetic origin are taken into the mechanics domain to calculate the deformation and vibrations, where the forces are fed as an input to the solid mechanics calculations. Until the computation of vibrations, the FEM is used, and for the acoustic noise, the BEM is employed. In the first step, the electromagnetic simulation gives the magnetic forces due to Maxwell stress, and these forces are then fed to the structural computation as an input excitation. The supply frequency was 1008 Hz, and three time periods were simulated in the electromagnetic 2D computation. There were 19,750 linear triangular elements in the 2D mesh, as shown in Figure 3a, and the structural 3D mesh of the stator contained 104,931 tetrahedral elements, which is given in Figure 3b. Electromagnetic and structural mechanics time-stepping simulations were performed with a time-step length of 2.5 e−<sup>6</sup> s for three time periods corresponding to a 0.003-s machine running time. The 2D electromagnetic computation took 17 min of CPU time to finish the simulation, and the 3D structural mechanics simulation took 4 h for each time period. A time-to-frequency FFT was done for the results of the third time period to transform the solution from a time to frequency domain. The acoustics computation to calculate the sound pressure level took 11 min of CPU time for each frequency. The magnetic flux density distribution across the cross-section of the motor is given in Figure 4.



**Figure 3.** The finite element mesh of the motor: (**a**) the 2D triangular mesh of the rotor and stator and (**b**) the 3D tetrahedral mesh of the stator.

**Figure 4.** Magnetic flux density distribution in the motor cross-section from the electromagnetic finite element method (FEM).

In the structural simulation, displacement and acceleration of the stator body were computed. The magnetic forces were computed using Equation (2), and these forces were fed as the force term into Equation (3) to calculate the displacements. As a second stage of structural analysis, a Fourier analysis of the results was performed to discover the major frequency components in the deformation spectrum. The deformation pattern of the stator is depicted in Figure 5, and the Fourier analysis results of the displacement at a point on the stator outer surface are given in Figure 6, where the rotor rotational frequency is 1000 Hz, the twice-supply frequency, 2*f* s is 2016 Hz, and 4032 Hz corresponds to 2*p*\*2*f* s, where *p* is the number of pole pairs.

**Figure 5.** Deformation of the stator due to magnetic forces from the solid-mechanics FEM. Displacement scale factor: 5000.

**Figure 6.** Frequency spectrum of the radial displacement of a point on the stator outer boundary.

In the final step of acoustics, the accelerations from the FE mechanics domain are imposed at the stator boundaries, as is theoretically explained in Equation (16). This specific feature couples acoustics with the structural analysis for an acoustic–structure interaction. For the acoustics computations, the air surrounding the motor is modeled as an infinite void, where no specific geometry or meshing is required. The BEM simulation calculates the acoustic pressure in the surrounding air of the motor and the resulting sound pressure level corresponding to different frequencies. Equation (15) is utilized in the calculation of sound pressure in the frequency domain. The acoustic pressure distribution outside the motor, corresponding to 2016 Hz of frequency, is shown in Figure 7. The sound pressure level produced by this acoustic pressure difference in the air is given in Figure 8. The same parameters for 4032 Hz are depicted in Figure 9a,b. Figure 10 shows the sound pressure level in far fields in different planes. The far-fields plots of the sound pressure level give a clear idea about the directivity of noise radiation in different planes around the motor at a particular frequency.

**Figure 7.** Acoustic pressure distribution in Pa outside the motor at 2016 Hz.

**Figure 8.** Sound pressure level in dB outside the motor at 2016 Hz.

**Figure 9.** (**a**) Acoustic pressure distribution in Pa outside the motor at 4032 Hz; (**b**) sound pressure level in dB outside the motor at 4032 Hz.

**Figure 10.** Far-field sound pressure level in dB at a distance of 1 m in the (**a**) *x*-*y* plane, (**b**) *x*-*z* plane, and (**c**) *y*-*z* plane.

#### **4. Discussion**

The results presented in this paper showcase how magnetic forces produce acoustic sound in an electrical motor. The emphasis of this study was mainly centered on a sound level computation through FEM–BEM coupling, and the results show that this multilevel approach combining two numerical methods and three different physics in a single finite element tool is an effective tactic for acoustic analysis of electrical motors. The structural mechanics results provide information about the deformation and vibration behavior of the machine and also the major frequency components in the vibration spectrum, as is shown in Figures 5 and 6. The computational results of the acoustic pressure variation occurring in the neighboring medium of the motor (Figures 7 and 9a) shed light on the pattern of sound pressure and how different frequencies cause dissimilar forms of pressure variations in terms of their directions and intensity, thereby varying the levels of sound (Figures 8 and 9b). The far-field sound pressure levels portrayed in Figure 10 offer an idea about how motor vibrations produce and direct sound in different directions.

One major drawback concerning numerical computations of acoustic noise is the relatively higher memory requirements and computational time required compared to analytical and semi-analytical models. The increased computational time of the FEM has been mentioned as a drawback of numerical methods in the literature, such as in References [7–9]. The use of the BEM in an acoustics domain fixes this issue to a considerable degree, without compromising the accuracy of the results. Although methods such as the permeance/mmf model [7,25] and the subdomain model [7,8,26] reduce the computational time, they possess some drawbacks in the modeling of complex geometries and in terms of the 3D effects of machines. In a multiphysical simulation environment, the use of the BEM instead of the FEM as the numerical method thus facilitates a competent framework for vibro-acoustic computations without compromising precision in detailed modeling, and it also has the benefit of reducing the computational time. If the FEM had been employed in the acoustic computations of the study presented in this paper, the entire surrounding air would have been meshed, and this would have caused increased memory requirements and computational time. A quantitative study comparing previous works could not be done because of the difference in motor types. However, qualitatively, the FEM–BEM approach has all of the benefits of numerical computations, especially in terms of accuracy and flexibility in modeling complex structures, and lessens the computational time significantly compared to a complete FEM model.

This paper does not include the measurement results of sound levels, and hence a future component of this study will be to conduct laboratory experiments to measure sound levels. In addition, to compare the measurement results to the simulations, the entire structure of the motor, including the frame and bearings, needs to be modeled in the future. Furthermore, computing and measuring the A-weighted sound levels corresponding to different frequencies need to be carried out in the second part of this study to precisely illustrate the sound pressure level in terms of human audibility [27]. In addition, in real-world problems where complex systems need to be simulated, materials cannot necessarily be connected: rather they are glued or clamped. Furthermore, the endplates and the frame of the motor could have an effect on the computed vibrations and sound pressures. Hence, the boundary layers in the structural FEM and acoustic BEM become more difficult to model. In those cases, accurate material models are required, and the model size could be increased. Both issues require measurements and iterative parameter adoption. However, the computational method presented in this paper facilitates a platform for vibro-acoustic studies by effectively modeling the acoustic noise of the motor.

#### **5. Conclusions**

The sound produced by an induction motor due to electromagnetic forces was successfully computed using a joint finite-element and boundary-element-based numerical computational technique. The energy conversion process in an electrical motor was modeled by showing how the electrical energy is converted into acoustic energy. The paper successfully illustrated how the proposed coupling methods and the FEM–BEM combination work effectively in acoustic studies of motors. The task of formulating a computational framework by combining theoretical equations in the right way and order, preparing numerical simulation stages, and coupling them to yield factual results in different stages and different physical domains was done in an original way in this study, as was presented in the paper. The proposed scheme is a vivid method describing how electromagnetic, mechanical, and acoustic domains can be analyzed and coupled in a numerical simulation platform. It provides an expedient and effectual tool for researchers in electrical machine and acoustics fields conducting magneto-mechanical and vibro-acoustics studies.

**Author Contributions:** Conceptualization, S.S.; methodology, S.S. and U.A.; software, S.S. and U.A.; validation, A.B., S.S., and U.A.; formal analysis, A.B., S.S., and U.A.; writing—original draft preparation, S.S.; writing—review and editing, S.S.; visualization, S.S.; supervision, A.B. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Conflicts of Interest:** The authors declare no conflicts of interest.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
