A Complete and Fast Analysis Procedure for Three-Phase Induction Motors Using Finite Element, Considering Skewing and Iron Losses

Featured Application: The aim of the work is to describe and propose a new approach for the induction machine analysis. The application of the method is immediate for fast and computational efﬁcient procedures for the design optimization and precise performance computation. Abstract: This paper deals with a complete ﬁnite-element analysis procedure for squirrel cage induction motors, including the presence of skewing and the iron losses evaluation. The machine is analyzed performing only magneto-static ﬁnite element analyses. Saturation phenomena are carefully considered in any operating condition, avoiding long time-stepping analyses. The synergy between analytical and ﬁnite element model leads to a rapid and precise estimation of the rotor induced current, saving computational time. Furthermore, the procedure proposed in this paper allows the motor performance to be directly derived, without the preliminary knowledge of the machine equivalent circuit. In order to complete the analysis, skewing effect is included, using the 2-D multi-slice technique, based on static simulations. Experimental tests are carried out and reported in order to verify analysis results.


Introduction
Despite the conventional structure of the Induction Motor (IM), the study of electromagnetic phenomena occurring in the machine is not immediate. For this reason, the most common way to obtain motor performance is using Time Domain (TD) Finite Element Analysis (FEA) [1]. The IM still plays an important role in the market, and accurate tools for its analysis and design are of great interest across industries [2,3].
To avoid a long computation time of time-stepping analysis, IM is often analyzed using the equivalent circuit [4,5]. The FEA can improve the parameters estimation, considering the skin effect in the rotor bars, and the saturation of the magnetizing inductance [6].
In this work, an accurate procedure is proposed: IM performance is derived by performing Magneto-Static (MS) FEA, in which both stator and rotor currents are imposed. In particular, the rotor current is computed connecting the FEA model with the inverse-Γ equivalent circuit, related to the Rotor Field-Oriented (RFO) model [7]. Basically, in the RFO analysis technique, stator and rotor d and q currents have to be imposed in order to verify the condition λ rq = 0, i.e., to reproduce the actual field distribution inside the machine and perform on-load MS simulations.
These features are included in the method described in [8]; in this paper some improvements to the MS IM analysis are introduced together with practical application examples.
In particular, the stray losses computation is not a trivial problem for the induction machine also supplied by a sinusoidal voltage source. In the classical FEA-augmented equivalent circuit, the iron losses resistance is derived by performing the no-load test in a single rotor position, taking the flux density peak value in the stator teeth and yoke.
Furthermore, as far as the stray on-load losses are concerned, the classical methods neglect this effect or use an equivalent resistance based on measurements. The proposed technique also includes the computation of iron losses due to the stator and rotor slot harmonics, which occur both during no-load and on-load operations, by means of the element by element technique [9,10]. Furthermore, the classical analysis does not consider the iron saturation due to the torque current, that can be serious during overload condition.
The skewing effect is considered. The skewed motor is represented by a number of 2-D slices, each at a different axial position in the machine. The quantities (e.g., torque and fluxes) are computed for each slice independently and then summed [11,12].
In [11], a 2-D multi-slice model is presented for the skewing analysis, avoiding the whole 3-D problem. Each 2-D FEA problem is solved in the time domain, imposing that each rotor slot has the same current in the considered slices, leading to a very large FEA problem.
In [12,13], an alternative procedure is proposed, combining the analytical voltage equation of the considered machine with FEAs, used to upload, in each time instant, the inductive lumped parameters to the analytical model, considering the skewing.
In this work, instead of TD analysis, MS problems are solved in each slice, since the rotor current has been already determined using the rotor field-oriented technique.

Rotor Equivalent Three-Phase Winding
To achieve a low harmonic content in the Magneto Motive Force (MMF), the stator winding of a common IM is well distributed along the periphery with a sufficient number of slots per pole per phase.
In such a case, only the first harmonic of air-gap flux can be considered, and spatial behavior of induced voltage and current, within the rotor bars, exhibits a sinusoidal waveform, as shown in Figure 1. . Induced voltage and current spatial distribution in the slots of a squirrel cage rotor when the air-gap flux density B g is sinusoidal. v sl indicates the relative speed between the stator field and the rotor; z bar = r bar + jx bar is the bar impedance; i bar and e bar are the bar-induced voltage and current; L stk is the stack length.
A useful trick is considering a three-phase rotor equivalent winding sinusoidally distributed in the rotor slots, in order to properly set the rotor current and compute the rotor flux linkages. The complete procedure to design such a winding is described in [8]. For low rotor frequencies, the rotor bar current distribution can be considered uniformly distributed.
The number of conductors per phase of the rotor winding is fixed to have the same number of effective conductors as the stator winding: This equivalence is made to facilitate the rotor current and flux linkage computation, since no transformation coefficient is needed. The condition (1) makes the stator and rotor magnetizing component of the synchronous inductance to be equal, since the two windings share the same magnetic path for the magnetizing flux: where g is the magnetic air-gap, which considers saturation and the Carter coefficient. Using this three-phase equivalent winding, it is possible to apply the Park transformation even for cage winding deriving d and q quantities in the RFO reference frame.

IM Inverse-Γ Equivalent Circuit and Its Connection to Finite Element Model
The challenge of simulating the IM under load, using MS FEAs, consists of the rotor current prediction outside the FEA problem, since the time derivative of the magnetic vector potential is not included in the formulation.
The equivalent circuit is a common tool to estimate the stator and rotor current space vectors and the motor performance [4,5]. However, the inductances do not take into account the iron saturation due to the presence of stator and rotor torque currents.
Anyway, the analytical model gives useful information about the relative position and the amplitude of current space vectors, which can be properly imposed in static FEAs. The RFO model is based on the inverse-Γ form of the equivalent circuit in Figure 2. In this model, the rotor flux space vector (referred to the stator) λ rs is determined only by the magnetizing stator current i µ . It lies along the d-axis of the steady state space vector diagram (Figure 2b). The magnetizing current is also the d-axis stator current i sd in the RFO model.
L t and L ϕ in Figure 2a are, respectively, the global leakage and the magnetizing inductance in the inverse-Γ model: where L r is the rotor three-phase rotor equivalent winding synchronous inductance. Stator and rotor voltage equations, related to the RFO model are: where R s and R rs are the stator and rotor resistance; ω s is the stator angular frequency. In the RFO model, the rotor current i rs represents the torque current i τ , which lies along the q-axis of the space vector diagram reference frame, Figure 2b. Introducing the transformation coefficient (L m /L r ), stator and rotor currents are: The rotor current space vector lies along the q-axis of the RFO reference frame. Using (5), stator and rotor flux linkage space vectors can be expressed as: where, in the RFO synchronous reference frame: The flux linkage-current relationships (7) show the first advantage of this IM model: the current i sd is involved in the main flux linkage creation, whilst the current i sq produces leakage flux linkage and implicates the torque production. Furthermore, rearranging (5), it is: The relationship (8) represents another peculiar feature of this model: the rotor current vector acts only along the q-axis of the RFO reference frame. This condition represents an advantage to create a static finite-element model of the machine. As shown in Figure 3, the dq reference frame is geometrically placed in an arbitrary position, which corresponds to the phase A magnetic axis for simplicity. In this reference frame, the stator current is imposed. From the analytical model, the rotor current has to be set along the q-axis and the amplitude has to follow the rule (8).
The condition (8) is related to the model feature λ rq = 0. So, the proper rotor current amplitude is set when the rotor q-axis flux linkage is equal to zero. The procedure to compute the proper rotor current, using MS FEAs will be described in detail. The final result of this approach is an instantaneous electromagnetic on-load working condition of the motor.
With the proper rotor current set in the MS problem, together with the stator current, the machine flux linkages can be obtained from the field solution. The mean torque is computed as: Furthermore, starting from the rotor voltage equation, the slip is derived as: Figure 3. In MS FEA, stator and rotor currents have to be imposed according to the inverse-Γ model: the condition λ rq = 0 has to be verified in order to properly analyze the machine (the rotor flux linkage space vector lies along the d-axis).

FE Analysis Procedure
The method described in this paper can be implemented using any FEA software. The MS FEA problems do not require particular attention in the mesh definition, since no induction phenomenon has to be observed. Normal triangular mesh is used with first order elements. To get more reliable results, the machine air-gap has to have at least three elements in the radial direction, and the iron regions, near the air-gap, to have a higher element density.
The fundamental condition to find the proper rotor current to set in the MS FEA model is defined in (8). If the rotor current follows this relationship, the flux λ rq is automatically set to zero, which is the checking result to verify that the rotor current amplitude, set in the simulation, is correct.
In MS FEAs, the rotor current is imposed using the equivalent rotor winding, introduced in Section 2. This is a three-phase winding, then similar to the stator, but the conductor allocation in each slot leads to a sinusoidal distribution of the current along the rotor periphery.
Considering again condition (8), the rotor current depends upon the ratio (L m /L r ), while i sq is imposed. The following FEA procedure allows a precise estimation of the inductance, carefully considering the iron saturation.

Step 1: Machine Parameters Estimation
Let us impose the stator d and q-axis currents and assume the dq reference frame position as in Figure 3. The rotor current, according to the inverse-Γ model, is set along the q-axis, and its amplitude is given by relationship (8).
The inductances L m and L r are unknown at the beginning: a first MS FEA is necessary to derive their value, making an initial hypothesis about i rq amplitude.
Since L m L r , a reasonable approximation is: i rq = −i sq . On the other hand, the rotor d-axis current is already set to zero.
In the first simulation of the procedure, stator and rotor current vectors are: At this step, the conditions λ rq = 0 are not verified. However, such analysis is necessary to derive stator and rotor leakage inductances L σs and L σr , besides the mutual synchronous inductance L m . In fact, substituting (11) in the flux linkages equations: The first field solution is shown in Figure 4. In Figure 4a, it can be noticed that the flux lines within the rotor widely cross the d-axis, meaning that the q-axis rotor flux linkage is not zero ( λ r is out of phase with respect to the d-axis). The flux linkages in (12) are related to the equivalent circuit in Figure 4b, where the mutual coupling is modelled with the T-form inductance network. The number of effective turns, chosen for the rotor equivalent winding, is the same as the stator. Thus, the transformation ratio n T = (N s k ws )/(N r k ws ) = 1.
From the field solution reported in Figure 4a, stator and rotor flux linkages are derived by integrating the magnetic vector potential within the slots. The machine inductance are derived as: These parameters are effective only in one specific working point, since, varying the motor current, or the rotor position, also the saturation map changes.

Step 2: Rotor Current Computation
The first step of the procedure ends with the computation of the T model inductances and, applying (8), the proper rotor current is achieved, according to the RFO model.
A second MS simulation is performed in which the stator current vector is the same as the previous step: only the rotor current is changed: Furthermore, in the second step, the rotor current space vector is imposed along the q-axis, but the amplitude has been modified.
The second simulation field solution is shown in Figure 5a. The spatial distribution of the current is the same as the previous simulation in Figure 4a. However, the amplitude of the rotor current is imposed according to the RFO model. In Figure 5a, the rotor flux lines are almost parallel to the d-axis, meaning that the q-axis component of the rotor flux linkage is close to zero. In Table 1, the procedure is applied on one-pole pairs, 11 kW induction motor, at different working points.  The equivalent circuit, related to the second MS FEA, is shown in Figure 5b. When the rotor current is inherent with the RFO model, the overall machine leakage flux is related to a single inductance, achieved as: whereas the magnetizing inductance L ϕ can be derived using the second in (3). From this second simulation, the motor torque and slip can be derived using (9) and (10).
Typically, in IMs, the ratio in (8) is almost in unity, since the leakage component of L r is a little percentage of the magnetizing component L m .
Actually, the rotor current imposed in the first step is a close approximation of the correct current according to the RFO model, imposed in the second step of the procedure; see Table 1, in which i sq = i rq corresponds to the rotor current in the first step and i rq is the rotor current in the second step.
For this reason, the parameters L m and L r computed in the first FEA, do not change significantly in the second FEA; in fact, the rotor current, computed by (14) and imposed in the second FEA, is already correct to set to zero the q-axis rotor flux linkage; see the comparison between λ rq in the two steps, reported in Table 1.

Skewing Analysis
The skewing effect analysis is performed connecting the multi-slice theory [11,12] with the procedure described above.
The continuous skewing is approximated as a stepped skewing, considering that in each slice the rotor is unskewed but properly rotated with respect to central reference position. For the induction motor also, the electric loading has to shift together with the rotor. The main meaning is that, with the skewing, the same bar current is distributed in a larger portion of space, depending to the skewing angle, interacting in different ways with the stator field in the different slices. The novelty here is that the rotor-induced current is predicted using the MS FEAs procedure, avoiding the time-stepping analysis.
In the following example, five 2D slices have been considered [14]. In this analysis, the IM exhibits a skewing angle of o slot. The angle of each slice, with respect to the dq reference frame, is defined as: where α sk is the rotor skewing angle, M is the considered slice, M tot is the total number of slices. The first considered slice is already reported in Figure 4a. In the other slices, the dq reference frame remains the same, together with the stator current vector. The rotor and its electric load are rotated backward or forward according to the slice angle α sl .
An example of a 2D slice is reported in Figure 6. The rotor skewing is represented by the reference frame α r β r rotation, forward and backward.
The field solutions are post-processed in the same way as for the slice aligned to the dq reference frame (α sl = 0): in particular, the angle between α r β r and dq, used in the Park transformation, to achieve stator and rotor λ d and λ q , is the same for each 2D FE problem. In other words, the rotor flux λ rβ r is considered valid even in the dq reference frame (λ rβ r = λ rq ) for each slice, and i rq = i rβ r . Basically, with skewing, the RFO Analysis method works in the same way as described in Section 4. In the first step of the procedure, stator and rotor current space vectors are imposed as in (11). A series of simulations is carried out, rotating backward and forward the rotor and its electric load.
Stator and rotor dq flux linkages are mediated to achieve values, taking into account the skewing. A first set of inductances can be derived from the field solution, using (13) and the average value of the fluxes.
Then, for the second step, the rotor current is imposed according to (14). A further set of slice simulations provides the new value of the rotor q-axis flux linkage; with skewing, the average value, considering each slice, has to be close to zero.
In Table 2, the analysis results are reported in several working points, increasing the torque. It can be noticed that, in the second step of the RFOA analysis, the rotor q-axis flux is significantly reduced. In Figure 7, stator and rotor dq flux linkage values are reported for each considered slice, imposing i rq = −i sq . In Figure 7a, regarding the stator flux λ sd , moving the rotor anti-clockwise, the rotor current partially acts along the d-axis, increasing the magnetizing flux. q-axis stator flux behavior is reported in Figure 7b: moving the rotor electric load forward and backward, the q-axis magneto-motive force balance is less effective and λ sq can increase.  As already mentioned, the rotor flux linkage λ rβ r in several slices is considered a part of the overall flux linkage λ rq . In other words, the transformation angle from rotor abc to dq is the same for each slice and it is defined by the angle between the rotor abc and dq in the central slice.
The flux linkage λ rβ r is linked with a theoretical winding with axis along β r . When a central slice is considered, this winding is able to gather even the d-axis flux, as shown in Figure 8.   In Figure 8b, a forward rotation of the rotor and the β r winding is considered, which corresponds to α sl < 0. The angle between the axes β r and d is called α, with In the case of α sl = 0, and i rq = −i sq , the flux λ rβ r is: where i rβ r = i rq and λ rβ r results negative, with L σr as the rotor leakage inductance.
With α < π/2 (Figure 8a), cos α > 0 and the contribution of Φ d to λ rβ r is positive: The rotor β r flux linkage can be written as: considering (17), the sum of the first two terms of (19) is negative, whereas the contribution of the d-axis flux is positive. Thus, for negative slice angles, the rotor β r -axis flux amplitude is decreasing, as can be observed in Figure 7c.
In Figure 8b, a backward rotation of the rotor is shown. In this case, the slice angle α sl is positive and α > π/2. As a consequence, applying (18), the cosine is negative and the d-axis flux contribution is negative: As before, the sum of first two terms is negative, as the contribution of Φ d . Thus, for positive slice angles, the flux λ rβ r is increasing its amplitude (see Figure 7c).
After the first set of simulations with i rq = −i sq , stator and rotor inductances are computed as: The rotor synchronous inductance is computed as: L r = L m + L σr and the rotor current for the second step of RFOA procedure is derived using (14).
The second set of slice simulations is carried out in the same way as before, imposing the uploaded value of i rq ; the rotor flux λ rq is plotted, for each slice, in Figure 9. . Rotor flux linked with the β r winding (λ rq = λ rβ r for each slice) for each considered slice, with the rotor current imposed according to (14).
From Figure 9, the mean value of λ rq is almost equal to zero. Furthermore, it can be ascertained that the rotor β r flux linkage variation is not symmetric with respect the central slice. With rotor skewing, it is like the dq plane position is not placed in the center of the stack length.
Finally, from the second step, the transient and magnetizing inductances, in the inverse-Γ circuit, are computed as: where L m and L r are achieved in the first step of the procedure.

Method Application 1: Average Torque and Current vs. Slip
The analysis procedure has been adopted for the simulation of a three-phase one-pole pair IM with cast aluminum rotor, nominal power P n = 11 kW, frequency f n = 50 Hz, phase voltage V n = 400 V. The motor has 30 rotor slots and 36 stator slots and a skewing angle of one slot pitch. The motor's efficiency class is IE3.
In order to achieve a good estimation of the motor torque, a complete simulation varying the rotor position is required, and, according to the position, the rotor current has to be recomputed to keep to zero the rotor q-axis flux, whereas the d-axis rotor current variation is neglected.
In this way, the torque ripple can be evaluated as well as the average value, imposing a certain slip and stator current: the d-axis current is fixed according to the rated voltage, whereas the q-axis torque current depends upon the slip. At the rated power, this motor works with s = 0.03, at 2910 rpm.
The torque behavior with and without skewing is shown in Figure 10. After the rotor skewing, the average value of the torque is reduced by 2.5% (it can be noticed that the T dq passes from 40 N m to 39 N m in Figure 10a,b), and the ripple passes from the 12.5% to 2.5% (it is evident looking at the behaviors of the Maxwell stress tensor torque with the rotor position).  In Table 3, the comparison of the torque harmonics in the unskewed and skewed motor is reported. The 6 th harmonic in the torque is due to the pulsation of the magneto-motive force amplitude, the other main harmonics are the stator and rotor slot harmonics and the combination of the two (high frequency ripple).
After the skewing of one slot pitch, the 30th harmonic is dramatically reduced together with the 180th, that is the l cm(36,30). The stator slot harmonic becomes one-half of the unskewed value.

Measurements Comparison
In order to validate the adopted IM model, measurements comparison has been carried out. The motor average torque and stator current have been measured in several working points, increasing the load from 25% up to 150% of the nominal value. The comparison with simulations is reported in Figure 11. A good agreement has been found. Figure 11. Measurements comparison about average torque and stator current. The measurements have been carried out according to the standard for grid-connected IMs' efficiency computation: keeping constant the voltage, the load is increased from 25 % to 150 % of the rated one.
Even if the effectiveness of this analysis has been proved, the number of required simulations is quite high, considering the high order harmonic in the torque behavior and the presence of the rotor skewing. Referring to Figure 10, a good estimation of the average torque comes from (9). The authors propose to compute the average torque of the motor using the T dq value, derived in a single rotor position, in order to save computational time.

Method Application 2: Iron Losses Computation
The IM analysis rotating the rotor allows for the variation of the flux density in each element, with the rotor position, to be derived. The flux-density behavior in the iron elements is then processed to derive the hysteresis and eddy current losses.
In Figure 12, the flux-density behavior in several parts of the machine is shown, and the presence of the high frequency component is evident. The fast variations of the flux density cause serious losses located especially in the stator and rotor tooth tips, as shown in Figure 13.  The high frequency components in the flux density are mainly due to the stator and rotor MMF slot harmonics and to the slot openings, in particular the stator ones.
There are two main methods for the losses computation with the element-by-element technique:

1.
Steinmetz method: Once the flux-density behavior during an electric period has been computed, the Fourier analysis is performed, and the Steinmetz formulation of the losses is applied to each flux density harmonic in each element [10].

2.
Bertotti method: Here, the hysteresis losses are imposed to be linked only to the main harmonic of flux density; on the other hand, the eddy current losses are computed considering the square of the flux-density time derivative. Such a method provides a sort of instantaneous value of the iron losses [9].

Steinmetz Formulation for Iron Losses Computation
The Steinmetz formulation has the following expression: with: γ lam lamination density V el element volume ν harmonic order k hy,ν hysteresis losses coefficient k ec,ν eddy current losses coefficient f s stator frequencŷ B el,ν ν-th harmonic flux density amplitude in the element The formulation is valid in the presence of both time and spatial harmonics; in this work, only space harmonics are considered.
Once the flux density behavior during an electric period has been computed, the Fourier analysis is carried out and (23) is applied to each flux density harmonic in each element.
The main advantage of this method is the harmonic division of the iron losses. However, at least one-half of the electric period has to be simulated to obtain the complete flux density waveform in each element of the mesh. The mirroring of one-half period waveform is made to obtain the one period behavior of the flux density required for the Fourier analysis.
Furthermore, due to the slip, the flux density in the rotor elements has a non-whole number of periods and the Fourier analysis is not correct. To overcome this problem in using this method for the IM, the authors propose neglecting the slip, imposing the rotorload current anyway. In post-processing, the rotor harmonics frequency will be correct considering the slip related to the first harmonic of flux: The different frequencies, induced by each spatial harmonic in the rotor, is due to the different numbers of poles that they have and the different speeds with respect to the stationary reference frame. This represents an approximation. Concerning the stator space harmonics, the correction in (24) is not proper: imposing the synchronism, e.g., the 5th and the 7th will have the same frequency, in the rotor reference frame, when actually they would induce different frequencies, when the slip occurs. This is true for all the stator magneto-motive force harmonics, including the slotting. In any case, this approximation is quite good if the slip is low as in the normal working condition of the IM.
The key point of this method is the choice of the hysteresis and eddy current losses coefficients. It is not trivial because flux density harmonics, with very different frequencies, coexist in the motor.
In particular in the IM, the fundamental stator harmonic and the stator and rotor slot harmonics are the main source of losses. The fundamental in the stator is 50 Hz; the first rotor slot harmonics in the stator have the frequencies: Thus, with flux density harmonics with such far values of frequency, it is not proper to consider only one couple of c hy , c ec coefficients for all the harmonics. The correct approach should be assigning to each harmonic a coherent couple of coefficient [15], interpolating the electrical steel datasheet with the closest frequencies to the considered one. The motor steel is a M400-50A, and the available specific iron losses datasheet includes the frequencies 50, 100, 200, 400, 1000 and 2500 Hz.
For the fundamental, the hysteresis and eddy current coefficients are derived, interpolating the 50 and 100 Hz losses curves at 1 T. Using these coefficients, the approximation of the datasheet curve at 50 Hz is shown in Figure 14a. The adopted Steinmetz-specific losses formulation is: Figure 14. Hysteresis and eddy current coefficient derived for the fundamental and the high frequency harmonics. The approximation of the datasheet-specific losses behavior is reported. (a) Hysteresis and eddy current coefficient for the fundamental flux density harmonic; (b) hysteresis and eddy current coefficient for the high-frequency flux density harmonics.
To derive the coefficients for the high frequency harmonics, the specific losses curves at 1000 and 2500 Hz have been interpolated at 1 T. The approximations of the datasheet curves, using the high frequencies coefficients, are reported in Figure 14b.
It can be noticed that the approximations in Figure 14b are good at least up to 1 T, and a greater amplitude of the flux density variation, due the slot harmonics, is not expected in this motor, except for a current much greater than the rated one.
In the IM, also the stator belt harmonics (5th, 7th, 11th, 13th, etc.) produce rotor iron losses, rotating, with respect to the rotor, at a higher frequency than the stator one, function of the harmonic order and the slip [16]. However, the belt harmonics' amplitude is remarkably lower than the slot harmonics and the frequency of the most important of them (5th and 7th) is not much higher than the fundamental. As will be shown later, the iron losses due to the stator belt harmonics are quite low and, for the sake of simplicity, the fundamental hysteresis and eddy current coefficients are assigned to these harmonics.
In Table 4, the stator and rotor harmonic losses are reported in the rated no-load and on-load operating points. Considering, at first, the no-load regime, it can be noticed that a serious amount of the total iron losses is due to the stator slot harmonics considering both the magneto-motive force harmonics and the slot openings. Furthermore, at no load, the stator losses due to the rotor slot harmonics, are almost negligible.
Moving to the on-load operating point, the harmonics in the rotor become non-integer, due to the presence of the slip. The iron losses due to the stator winding belt harmonics are linked to the harmonics 5.82 and 11.65, applying (24) (at no-load, the losses due to these harmonics, are negligible). For the same reason, the stator slot harmonics in the rotor become of order 34.96, 69.92, etc.
With the load current flowing in the rotor cage, the stator iron losses due to the rotor slot harmonics become an important amount of the total losses. Considering the current Standard definition, these kinds of losses are part of the on-load additional losses.
During the on-load operation, the first harmonic losses are lower than in the on-load condition because of the voltage drop on the stator impedance, considering that the stator is supplied by a constant voltage source.

Bertotti Formulation
The other main formulation for the iron losses, linked to the element by element analysis, is the Bertotti formulation: with: B el = B 2 x,el + B 2 y,el , instantaneous element flux density d lamination thickness σ lamination conductivity For this lamination, the excess losses coefficient, at the rated frequency, is equal to zero.
The instantaneous values of the stator and rotor iron losses are reported in Figure 15. It can be noticed that, using this method, it is not necessary to simulate a complete electric period; actually, simulating only one stator slot pitch is sufficient to get the average value. Furthermore, the comparison between the two considered methods is reported in Table 5.  Finally, the main advantage of the instantaneous losses method is the computation time: a much lower time span has to be simulated to get the average losses. The drawback is that the harmonic subdivision of the losses is not possible anymore. Another advantage is that the presence of the slip is not an issue, since no harmonic analysis is required.
The skewing analysis in Section 5 can be coupled with iron losses computation, to get the variation of the specific losses with the slice angle. In Figure 16, it can be verified that in the skewed machine the iron losses are almost equal to the unskewed case.

Measurements Comparison
At first, the comparison between measurements and simulation results concerns only the motor considered in Section 6. The no-load test on the machine has been performed applying the stator voltages 30%, 40%, 50%, 70%, 85%, 100% and 115% of the rated value. Then, the iron losses are computed, following the procedure reported in the Standard for the AC machines' efficiency evaluation. The first comparison result is reported in Figure 17. Even if the model takes carefully into account the role of the high frequency harmonics in the losses computation, the simulations do not fill very well the measurements. Possible reasons are related to the implemented losses model: in (23),B el,ν represents the amplitude of the flux density harmonic, thus the model does not consider the presence of rotational fields [17,18], that occurs specially in the top of the stator teeth and in the stator and rotor tooth tips [19]. A further improvement can be obtained by summing the losses components, related to B x,ν and B y,ν , in each element.
Another possible reason is that the adopted model does not take into consideration the manufacturing defects introduced by the punching and the rotor turning.
Generally, the increasing factor is applied to the overall value of the calculated losses, including also the first harmonic component. Actually, as shown in Figure 13, the high frequency losses are placed where the manufacturing defects can have a strong influence: in the stator and rotor teeth tips, where, in addition, an high concentration of rotational field losses take place [19]. On the other hand, the first harmonic flux lines are parallel (in the stator teeth and back iron) to those surfaces that have been punched. For these reasons, it can be defined an increasing factor that concerns only the high frequency components of the losses. Considering the 11 kW prototype, for each voltage level, an increasing factor for the iron losses has been defined as: where P meas are the measured iron losses, P sim,1 are the first harmonic losses from the simulation and ∑ ν>1 P sim,ν are the high frequency losses from the simulation. It has been found out that the value of k i is almost constant for each voltage level, except the last one, where it is rather higher.
For the highest value of the voltage, the flux density in the iron reaches the saturation level (about 1.75-1.85 T) and it is worth to notice that the approximation in Figure 14a was able to fit neither the datasheet curve at high flux density, using the formulation in (27). Thus, the increasing factor is defined considering the voltage levels lower or equal to the rated value (400 V); making the average value, k i = 2.2.
Furthermore, the same increasing factor for the high frequency losses has been used for several motors in order to test the validity for different prototypes.
Several motors belonging to the frame size IEC 132 have been measured, with stack length 150 mm, 180 mm and 200 mm, with a mechanical power of 7.5 kW, 9.5 kW and 11 kW, respectively. For all the motors the lamination is the M400-50A. After the correction of the high frequency losses, the comparison with no-load measurements for the several motors is shown in Figure 18. During the on-load operation, the stator iron losses increase due to the rotor slot harmonics, produced by the rotor magnetic load. Furthermore, in the cage, the induced current due to the stator harmonic fields produces additional Joule losses. The on-load additional iron losses are computed as: where P fe,OL are the on-load iron losses and P fe,NL are the no-load iron losses. The additional cage losses can be estimated analytically, using the harmonic equivalent circuits. From measurements, the additional losses are computed according to the Standard. Comparison with the simulations results is reported in Figure 19. Finally, the skewing influence on iron losses can be neglected, as suggested in [20].

Discussion and Conclusions
The method described in this paper allows a very fast induction motor analysis using magneto-static finite element analyses. The rotor current is computed linking the analytical laws in the rotor field-oriented reference frame with the precise estimation of the parameters from the finite element model. The result is a rapid procedure which avoids the time-domain formulation to get a precise estimation of the performance.
The time-domain analysis, run with any commercial software, requires many minutes or even hours to get the steady state solution, especially because of the electrical transient. Furthermore, in the presence of rotor skewing, the finite element problem becomes more complicated and the computational effort much greater. The solution proposed in this paper adopts a hybrid analytical finite-element approach which allows the immediate computation of the rotor current and the easy implementation of the skewing effect. Thus, the computation time to achieve steady state performance requires less than a minute with the proper mesh and when applying the correct periodicity. The classic analytical model results in being the faster alternative, but the saturation model is too simplified and does not consider the contribution of the on-load current. Furthermore, only simple iron losses models can be coupled with the analytical equivalent circuit analysis.
The method is useful also for skewed machines, representing a computational efficient approach. The final comparison between experimental and simulation results proves the effectiveness of the new procedure.
The applicability of this analysis method is actually rather broad. Such a simple and rapid analysis method is very useful in the machine design stage. In fact, to the aim of achieving the optimal motor geometry, a huge number of solutions could be analysed in a relatively low time.
The first example application reported in this work concerns the computation of the torque and current vs. slip characteristics, obtaining good results compared to experiments.
The second application method shows the possibility of computing the iron losses whilst accounting for the contribution of the stator and rotor space harmonics. The procedure is based on the storing of the flux density behavior in each element, considering a central rotation. The losses are computed in post-processing, analyzing the flux-density variation.
The proposed method can find application as the fast analysis procedure in a design optimization tool. The proposed approach can be useful, especially in the automotive sector, where special induction motors have to be designed, not following the "industrial design rules".
In conclusion, the induction motor analysis procedure proposed in this paper represents a valid alternative to the common time-domain analysis. The method is actually based on the equivalent circuit theory, but the computation of the flux linkages carefully considers the iron saturation in every working condition.