Transient Modeling of Induction Motors considering Space Harmonics

This paper presents a transient mathematical model for induction motors based on space harmonics and its validation by finite-elements method comparison. The model is based on the computation of the magnetomotive force-wave for squirrel cage induction motors, and accounts for the geometry and physical dimensions of the stator and rotor windings. The harmonic components of the airgap field are calculated at each time-step considering the currents in the stator windings and in the rotor bars, the latter modeled as overlap of contiguous coils, as well as the rotor position. The proposed model is compared with finite elements calculations to outline differences and effects of simplifications and harmonic order. A 200kW induction motor designed for a premium electric vehicle is assumed as case-study.

e.m.f. induced in the ݄-th stator winding; ݈ ఙ௦ , ‫ݎ‬ ௦ leakage inductance, resistance of the stator winding; ݈ ǡ ‫ݎ‬ rotor bar inductance, resistance; ݈ ǡ ‫ݎ‬ rotor end-ring inductance, resistance between consecutive bars; ‫‬ pole-pairs; ܰ ௦ number of series conductors per slot; ‫ݍ‬ ௦ stator slots per pole and phase; ‫ܬ‬ load and rotor inertia; ‫ܮ‬ equivalent length of the machine; ܰ number of rotor bars; ܴ mean radius of the air-gap; ‫,݉ܫ‬ ܴ݁ imaginary and real part; ߜ ఔ equivalent air-gap related to the ߥ-th harmonic; ߴ, ߱ electric rotor position, speed; ߶ ത ఔ ߥ-th order space vector of spatial distribution ߶; ߤ air permeability; ‫ܭ‬ ఔ winding factor related to the ߥ-th harmonic; In recent years there has been a renewed interest in the dynamic modeling of Induction Motors (IM), also due to their use in new and particularly relevant application areas, such as automotive [1][2] [3] .
Lumped Parameters (LP) models are the most used method to evaluate the transient performance of electrical motor drives, [4]. Based on a set of differential non-linear equations with constant or variable parameters, they allow a dynamic analysis with small calculation times with obvious advantages. Unfortunately, the classical analysis that uses LP models, does not take into account the air-gap space harmonics.
Indeed, the space harmonics induce voltages and harmonic currents circulating in the rotor windings. The interaction between the harmonic currents generated in the rotor and the harmonic fluxes results in harmonic torques, vibrations and the noise. These are some of the reasons of the so-called stray torques and additional losses, and in some cases they are sufficient to cause the speed-torque cusps that result in subsynchronous crawling [5].
Space harmonics are mainly considered by motor designer [6] but in recent years they are becoming more familiar to power systems engineers dealing with modern drive, due to the increase of variable frequency power systems applications [7]. Different methods have been developed in literature to take into account the space harmonics in transient analysis [8][9] [10].
The main characteristic of these methods is the need to set the LPs of the model considering an off-line FEM model, with the evident disadvantages of a non-analytical implementation in terms of computational load and detuning of the FEM model.
The use of co-simulation tools based on finite elements method (FEM) is recognized as the most precise way, until today, for modeling the coupled field circuits and motion of IMs, accounting for both saturation effects and space harmonics, [11]. Indeed, the modeling of the mechanical motion and field sources variation simultaneously allow coupling the instantaneous stator and rotor fields in transient behavior.
The main drawback of co-simulation is the execution time consumption due to the necessity to frequently recall the computationally heavy FEM model to guarantee an appropriate ratio between the calculation time-step and element size and to remesh the air-gap [12].

Transient Modeling of Induction Motors considering Space Harmonics
This paper presents a transient model suitable for the simulation of cage induction motors accounting for the actual topology of the machine, based on the space vector theory. The model adopts an analytical approach to evaluate the Fourier components of the airgap Magneto-Motive Force (MMF), called as MMF-Wave (MMFW) in the paper.
The equations account for the geometry and physical dimensions of the stator and rotor windings and the rotor position considering the field sources, represented by the currents in the stator windings and in the rotor bars. The analytical develop is presented in section II.
In section III, the MMFW approach is compared with FEM computations in order to evaluate the impact of simplifying hypotheses on the overall accuracy. The case study is a 200 kW induction motor designed for traction of an electric vehicle in the frame of the Horizon 2020 "Rare Earth Free Drive" project.
Finally, in section IV, the transient model implementation is presented, based on the numerical integration of the stator and rotor voltage equations calling the MMFW procedure in every calculation time-step. Comparative elements are given in terms of calculation time with respect to the transient FEM method.

III. DYNAMIC MODEL
The proposed IM model is developed on the assumptions of uniform and radial air-gap magnetic field and infinite iron permeability.
The machine is represented in terms of the three-phase equivalent circuits for the stator winding and polyphase circuits for the rotor one. By assuming the stator and rotor currents as state variables, the differential model is given by by the voltage balances as follows: where ܼ ௦ ൌ ሺ‫ݎ‬ ௦ ݈ ఙ௦ ௗ ௗ௧ ሻሾ‫ܫ‬ሿ ሺଷ௫ଷሻ and ܼ is an equivalent full matrix for the rotor which accounts for the resistance and leakage inductance of the rotor bar ‫ݎ‬ and ݈ and respective section of the cage end-ring ‫ݎ‬ and ݈ (see Fig. 3 and its explanation in the following of this section).
Assuming that internal structure of the winding is symmetric only space harmonics by odd order are present. The elements of the e.m.f.s. matrices ሾ݁ ௦ ሿ and ሾ݁ ሿ according to the space vector theory can be expressed as: ܰ is a coefficient introduced to evaluate the flux linkage by each rotor bar (see Appendix), and ߶ ത ଶఔାଵ is the (ʹߥͳ)-th order space vector of yoke resultant flux given by: with ߪ ത ௦ǡଶఔାଵ and ߪ ത ǡଶఔାଵ the stator and rotor linear current density, respectively, and The mathematical model of the motor is completed by the mechanical equations: where the electromagnetic torque ܶ is given by: with ܶ that represents the load torque and ଓҧ ௦ defined in (9) By means of the numerical integration of (1), (2), (6) and (7), it is possible to calculate the stator currents and the electromagnetic torque.

A. Stator model
The stator flux in (5) can be calculated by the classical formulas known in literature for winding distributed in regularly spaced slots (symmetric arrangement) along the airgap, assuming an infinitesimal slot-opening width i.e. impulsive distribution of the linear current density.
For the case of study, the single-phase winding is realized with three identical full-pitch winding sections connected in series. The coils are positioned inside three contiguous slots as presented in Fig. 1. The following expressions are obtained for the linear current density and magnetic fields harmonic space vectors, respectively: B. Rotor model Due to its structure, the rotor winding of squirrel cage IM cannot be represented as the three-phase stator one if not in equivalent terms therefore losing the geometric detail.
In fact, at any time during operation, the current in each bar of the cage is different from the currents in the other bars as shown in Fig. 2. Moreover, the current flows in the end rings too. The matter is addressed in literature looking at the squirrel cage as an atypical case of polyphase winding, [13].
In particular, the squirrel cage can be considered as a polyphase "ring" winding with a phase number equal to the number of bars. It follows that the current in each bar can be computed as the difference of consecutive phase currents as shown in Fig. 3 and the induced voltage for each rotor mesh (phase) can be represented as a voltage source located in series with the end-ring components. The magnetic effects due to the proposed winding structure can be analyzed considering separately the magnetic field produced by a rotor phase current (i.e. the current flowing in opposite directions on a couple of consecutive bars).
By analogy of the three-phase stator winding and assuming the same hypothesis of infinitesimal slot-opening width, the MMF waveform due to each rotor phase can be represented as shown in Fig. 4: Fig. 4 Airgap MMF waveform produced by a rotor polyphase current As the figure shows, the flux produced by the elementary winding ߛ ൌ ʹɎȀ wide passes through the airgap and returns spanning the angle ɔ ൌ ʹɎ െ ɀ.
In order to evaluate the Fourier components ‫ܨܯܯ‬ ߥ of the resultant "rectangular" (with null mean value) MMF distribution in the airgap ݂݉݉ሺߙሻ, the well known normalized discrete ‫ܿ݊݅ݏ‬ሺ݂݊ሻ function can be borrowed, defined as ‫݊݅ݏ‬ሺߨ݂݊ሻȀߨ݂݊ with ݊ positive: The harmonic space vector of the linear current density, and magnetic fields, respectively becomes: where ߟ ߥ ʹߨ ܰ is related to the representation of the rotor MMF waveform through the ‫ܿ݊݅ݏ‬ function (see Appendix).

IV. CASE STUDY AND VALIDATION
The case study refers to the IM motor developed in the frame of the Horizon 2020 project ReFreeDrive "Rare earth Free e-Drives featuring low cost manufacturing" as the traction engine of high power 200kW electrical vehicles, [4].
Key points of such motor design were: 1) cost reduction through the minimization of the motor size; 2) optimized shape of the windings by hairpin technology, which permits to have a high efficiency (more than 94% considering WLTP class 3 drive cycle). Fig. 5 shows the cross section and main geometric data of the machine, while its rated point performance are summarized in TABLE I.  In this study the FEM model of the motor has been modified respect to the design set-up in order to match a few characteristics of the proposed space harmonic model and thus allow for a safe validation of this latter, and particularly: − infinite iron permeability is assumed; − rotor bar opening has been introduced to avoid a nonrealistic leakage inductance (see Fig. 6).

A. Validation analysis by static tests
That portion of the MMFW algorithm referred to field generation has been analyzed by static test, i.e. simulating the impression of constant current patterns in the motor windings. The comparison with the FEM model output in same conditions allowed both to validate the overall computation and evaluate the accuracy of the method. Results are presented in term of airgap magnetic field and include both distribution waveforms and Fourier spectrums.
Different supply conditions are evaluated, by feeding currents in the only stator windings, in the only rotor cage, or in both, as described in   The results show a good correspondence between the MMFW and the FEM computations, both with field sources in the stator (case A and B) and in the rotor (case C and D). The case E demonstrates the possibility to evaluate the magnetic field in the airgap by the superposition of the stator and rotor sources. The main differences between the two approaches are due to the hypothesis of assuming an infinitesimal slot-opening width. In fact, as well explained in [14], there are "field weakening peaks" in correspondence of both the stator and rotor slot-opening which are accounted by the FEM model but not by MMFW. This effect is clearly visible by the comparison of the magnetic field waveforms in the presented results and is usually managed by the Carter coefficient (and therefore defining an equivalent airgap as a function of the harmonic ߜ ఔ ). In this respect the mechanical airgap (0.8 mm) is purposely considered in the MMFW curves to allow for a reference evaluation of this effect.

B. Transient model implementation and tests
The dynamic model presented in section III, arranged for the electric equations as in (16), has been implemented in Matlab TM and solved by using a 4-th order Runge-Kutta (RK) integration algorithm, a good compromise between calculation time and accuracy of solution [15].
The block-scheme of the algorithm is shown in Fig. 12. The algorithm uses the information from one time-step to compute the next one ("single-step" algorithm). In each timestep the FFMW is recalled evaluating the differential inductances matrix ‫ܮ‬ ௗ and the vector that links the variation of the fluxes linkage respect the rotor position ௗథ ௗԂ . In particular, ‫ܮ‬ ௗ and ௗథ ௗԂ are calculated considering the incremental ratio of the linkage fluxes respect to a "small" variation respectively of the currents and the position, considering the currents and position values evaluated in the previous step. With the hypothesis of linearity of the materials, the ‫ܮ‬ ௗ matrix is constant. In the first test case the stator windings are excited by a voltage sudden step, with locked rotor. This test is performed to evaluate the transients of both the stator currents and the currents induced in the rotor bars avoiding disturbances due to rotor movement. Fig. 13 shows the stator currents transient and Fig. 14   As the figures show, the steady state values attained by the two methods are very close. This has been achieved thanks to the optimal approximation of the harmonic Carter coefficients used in MMFW following the observation pointed out in the previous section. The main differences can be observed in the transient behavior, affected mainly by variable rotor resistance due the skin effect introduced by the sudden voltage step.
The second test case refers to fixed rated speed operation (6000 rpm) with sinusoidal voltage supply and constant slip speed (337 ܸ ௫ stator voltage, 204 Hz frequency, 0.0339 slip speed). Fig. 15 represents the evaluated torque computed by the two compared methods in function of the time. In this case, the differences between the two modelling methods are more evident, both in the steady state and in transient dynamic. These can be partly ascribe to the skin effect due to the variable current in the rotor bars, nevertheless this aspect need further insight. The difference between the mean values of the torque at steady state is about 6,6 % and the rms values of the current of about 3 %.
Differently from the previous case the effect of the movement is simulated, hence the element mesh size of the FEM model and the time step must be reduced to have a good precision, with heavy effects on calculation times. To give an indication, the simulation time for the second test case is 3 second for the MMFW method respect to the 34 min for the transient FEM, considering a classical workstation (Platinum Intel® Xeon® 8253 2,20 GHz and 128 Gb of RAM).

V. CONCLUSIONS
In this paper, a transient mathematical model for induction motors based on space harmonics has been presented and validated by comparing finite-elements method.
The proposed analytical method makes use of simplifications to lay the theoretical foundations for the development of accurate and rapid transient analysis tools, reducing the huge computational cost of approaches based on co-simulation. The next developments concern the removal of the linearity hypothesis, considering the saturation of the materials, and adopting variable electrical parameters as, just to name a few, the bar and hairpin equivalent resistances.

VI. APPENDIX
Defined ‫ܫ‬ as the integral of the ݊-th harmonic with amplitude ‫ܣ‬ on a semiperiod: the integral ‫ܫ‬ ఊ of the same harmonic on a generic angle ߛ can be express in function of ‫ܫ‬ as follows: Since sine is an odd function, sin(a)-sin(-a)=2sin(a), Thereafter, the coefficient ߟ ఊ can be defined as the ratio between ‫ܫ‬ and ‫ܫ‬ ఊ :