Simulation of Selected Induction Motor Operating Conditions Using COMSOL Software

In this paper, some results from a simulation of three-phase squirrel-cage induction motor using COMSOL Multiphysics software are shown for the selected motor operating conditions. The simulation has been carried out under the nominal operating condition at the nominal rotational speed, for the operation under no-load and for the operation under locked rotor. The waveforms of the stator current, the currents in two neighbouring rotor bars, the current in the rotor ring and the motor inner torque are presented for each tested condition. The simulated values of some motor parameters are compared to their values from the motor catalogue, the manufacturer’s measurement, or the measurement in our laboratory, respectively. The simulated values of the rotor bar current and the rotor ring current are compared to their values from the scientific literature. The parameters of the motor equivalent circuit are calculated by COMSOL and compared to the values calculated from the motor data plate and from the measured values.


Introduction
To design, study, calculate or optimize the electrical machines, the Finite Element Method (FEM) can be used [1], [2] and [3].The FEM is a numerical technique which subdivides a solved domain into smaller parts called finite elements, which has several advantages [4].Each element is represented by a set of element equations.The solution is found in the single node points of each element.There are various types of software programs based on the FEM; here, the COMSOL Multiphysics has been used [5].
The induction motor was modelled as a twodimensional model with defined out of plane thickness, i.e. as the quasi-3D model.The Rotating Machinery, Magnetic interface under the AC/DC Module of COMSOL is used for its stationary and time-domain modelling [6], [7] and [8].This physical interface solves Maxwell's equations formulated using a combination of the magnetic vector potential and the magnetic scalar potential as the dependent variables.Under the Rotating Machinery, Magnetic interface the nodes Ampere's Law are used solving the magnetic field in the motor [9], [10] and [11].
The aim of this paper is a comparison of simulated motor parameters to its parameters from the motor catalogue and practical measurements, or to the theoretical presumptions, respectively.After this verification, the model will be used for studies of the tested motor under electrical and mechanical faults or anomalies.

Model of Induction Motor
The model has been carried out as a non-linear one with the 800-6AM type of electrical steel for the motor laminations.The B-H curve of this material was used in the Materials physics of COMSOL for both the stator and rotor laminations.Due to the laminations, homogenization approaches are applied for a determination of the anisotropic equivalent conductivity [12], [13] and [14].In this work, the homogenization approach was applied by the stated formulas: where σ is the conductivity of the iron sheet in (S•m −1 ), σ z is the equivalent conductivity in the direction perpendicular to the iron sheets, n is the number of the stacked sheets.
To feed induction motor from three-phase sinusoidal voltage source, the Electrical Circuit interface of COMSOL has been used, which includes connections to the Rotating Machinery, Magnetic interface and components to represent end winding impedances.This interface is also used for creating the external rotor circuit [8], [15] and [16].
The resistances of the stator winding overhang R Send , a rotor end-ring segment R Rend and a part of the rotor bar outside of the rotor iron R Bout are calculated using the stator winding and rotor cage parameters, and the cage geometry, respectively.The formulas are stated as: where ρ is the resistivity in (Ω•m), l 0 is the mean length of the stator winding overhang in (m), N S is the number of turns per one stator phase, A is the cross-section area of the stator winding conductor in (m 2 ), d RN is the mean diameter of the cage end-ring in (m), Q is the number of the rotor bars, A R is the cross-section area of the cage end-ring in (m 2 ), l Bout is the mean length of the bar outside of the iron sheets in (m), A B is the cross-section area of the bar in (m 2 ).
The temperature of the stator winding was detected from the motor measurements under nominal load and no-load, and taken into account for the resistivity calculation.
There are various approaches to calculate the motor leakage inductances as can be found, e.g., in [17], [18], [19] and [20].They depend, among others, on the determination of the leakage permeance factor.Whereas in [17], one common permeance factor of the stray around the winding overhangs is reported to be constant for almost all types of winding, in [18], it is stated separately for the stator and rotor.In [19], the endwinding leakage permeance factor is defined for various combinations of the stator and rotor winding types when calculating the total end-winding leakage inductance.Nevertheless, the formula for calculating the leakage inductance of the short-circuit ring of the cage winding is introduced in [19] as follows: where m S is the number of the stator phases, p is the number of the pole pairs, l B is the length of the rotor bar in (m), l S is the equivalent length of the stator in (m), v is the factor equal to 0.18 for machines with p > 1 as defined in [19].
In [20], the end-winding leakage permeance factors are calculated taking into account the various types of winding, winding parameters and geometry.The leakage permeance factor of the stator winding overhang is calculated from the formula: where k 0 is a factor defined in [20] depending on the type of the stator winding, y is the coil span in (m), q is the number of slots per pole and per phase.
The leakage permeance factor of the cage winding is stated as: where τ p is the pole pitch, g is a factor defined in [20] depending on the geometry of the stator and rotor winding overhang.
Even if there are some differences in the statements or calculations of the leakage permeance using these approaches, they are of no great significance, resulting in only minor impact on the results of the induction motor simulation.
The leakage inductance of the stator winding overhang was calculated from the formula in [18]: and the leakage inductance of the rotor end-ring segment was calculated using the formula in [18]: The values of both leakage inductances calculated from Eq. ( 9) and Eq.(10) were used in the Electrical Circuit interface of the motor model.

Simulation of Nominal Operating Condition
The induction motor with the catalogue parameters in Tab. 1 has been simulated.
The motor nominal speed from Tab. 1 was adjusted in the Rotating Machinery, Magnetic physics of COM-SOL.In Fig. 1, the induction motor geometry is shown  illustrating the magnetic flux density distribution and equipotential lines of the magnetic vector potential.
From the results of the motor simulation, the stator current in one phase, the currents in two neighbouring rotor bars, the current in the rotor ring, and the motor inner torque are shown in Fig. 2, Fig. 3, Fig. 4, and in Fig. 5, respectively.The harmonic analysis of the simulated currents has been carried out to detect the first harmonic to be used in comparing currents.The first harmonic of one cycle of the stator current is shown in Fig. 6.The mean value of the inner torque at the motor nominal speed was detected from its analysis, and it is shown in Fig. 7.
The comparison of the stator current and motor torque from the simulation to the ones from the motor catalogue and from the measurement in our laboratory is presented in Tab. 2.
The currents of the rotor bar and ring are compared to the theoretical results.The current of the rotor bar can be calculated from the formula, as found, e.g., in [19]: The current of the rotor ring is calculated using Kirchhoff's first law at each connection point of the bar and the ring from the formula, as found, e.g., in The first harmonics of one cycle of the neighbouring rotor bar currents and rotor ring current are shown in Fig. 8, and Fig. 9, respectively.
The comparison of the rotor bar and ring currents from the simulation to the ones calculated from Eq. ( 11), and Eq. ( 12), respectively, is presented in Tab. 3.
As seen in Fig. 8, the time period of the bar current is equal to 0.5 s, which corresponds to the frequency value of 2 Hz.The value of the motor nominal slip is 0.04 corresponding to the bar current frequency f r of 2 Hz according to the formula: where s is the slip, f is the frequency of the supply voltage (50 Hz).
The phase shift between the waves of the currents in the neighbouring bars from Fig. 8 detected from their harmonic analysis is equal to the value of 0.449 rad (i.e.25.730) being identical to the value calculated from the formula, as found, e.g., in [20]:

Simulation Under No-Load Condition
The simulation was carried out for the motor speed equal to 1,498.5 rpm, which is the value detected from the motor simulation for its no-loaded start-up [8].The waveform of the current in stator is shown in Fig. 10.The value of the stator current under no-load is equal to 5.5 A, as it was detected from its harmonic analysis.This value is higher in comparison to the value of 4.76 A from the manufacturer's measurement, but close to the value of 5.7 A measured in our laboratory.

Simulation of Motor Under Locked Rotor Condition
The simulation under the locked rotor condition was carried out for the motor speed equal to 0 rpm adjusted in the Rotating Machinery, Magnetic physics.In Fig. 11, the magnetic flux density and equipotential lines of the magnetic vector potential are shown at the nominal supply voltage.As seen, the flux in the rotor is pushed towards the motor air gap due to the high currents in the rotor cage.As a result, the flux is prevailingly accumulated in proximity to the rotor surface and between the slots, which are highly saturated.The magnetic flux density in the rest of the rotor iron is low, as it can be seen in comparison to Fig. 1 under the motor nominal operating condition.
The ambient laboratory temperature was taken into account for the simulation under the locked rotor condition.The stator current is compared to the catalogue starting current at the motor standstill when the motor temperature is assumed to be equal to the ambient one.The manufacturer's measurement was carried out using a quick-acting supply voltage under locked rotor   condition, so the motor temperature rise is not taken into consideration.
From the results of the motor simulation, the stator current in one phase, the currents in two neighbouring rotor bars, the current in the rotor ring, and the motor inner torque are shown in Fig. 12, Fig. 13, Fig. 14, and Fig. 15, respectively.
The RMS value of the stator current detected from its harmonic analysis is compared to the value from the manufacturer's short-circuit characteristic and to the value of the starting current in the motor catalogue in Tab. 4.  The data for the manufacturer's short-circuit characteristic were obtained from measurements up to 90.2 % of the motor nominal voltage, and from an extrapolation for higher values, i.e. also for the nominal voltage.
It was confirmed from the harmonic analysis of the currents that the rotor bar and ring currents frequency is equal to the frequency of the stator current, as it can also be detected from comparison of the waveforms in Fig. 12, Fig. 13, and Fig. 14 respectively.

Parameters Calculated from COMSOL Analysis
The equivalent circuit of the induction motor in the harmonic steady state is shown in Fig. 16, neglecting the iron loss.The stator phase resistance R 1 and the stator phase leakage inductance L 1σ are composed of two portions given by the formulas: The first portion of the formulas Eq. ( 15) and Eq. ( 16) includes the circuit model as defined by Eq. ( 3), and Eq. ( 9), respectively.The second portion includes the 2D field model.The parameters R Send and L Send substitute in 2D modelling the end winding parameters of the stator overhangs, and R Sf em and L σSf em are the parameters of the stator windings inserted in the slots modelled by COMSOL physics.The resistance R Sf em is calculated automatically by COM-SOL because the slots are defined as the multi-turn coil domains.Note that equation Eq. ( 16) for the rotor can be written in a similar way.The inductances L σSf em and L σRf em , exactly their sum, has been determined indirectly by calculating the corresponding stator leakage linkage flux Ψ 1SL under the locked rotor simulation condition.

Source
Stator current Motor catalogue 50.4A Manufacturer's value 54.83A Simulation 55.5 A

Parameters Calculated From COMSOL Analysis
The equivalent circuit of the induction motor in the harmonic steady state is shown in Fig. 16, neglecting the iron loss.The stator phase resistance R1 and the stator phase leakage inductance L1 are composed of two portions given by the formulas: R1=RSend+ RSfem (15) L1σ=LSend+ LσSfem.(16) The first portion of the formulas ( 15) and ( 16) includes the circuit model as defined by (3), and (9), respectively.The second portion includes the 2D field model.The parameters RSend and LSend substitute in 2D modelling the end winding parameters of the stator overhangs, and RSfem and LσSfem are the parameters of the stator windings inserted in the slots modelled by COMSOL physics.The resistance RSfem is calculated automatically by COMSOL because the slots are defined as the multi-turn coil domains.Note that equation (16) for the rotor can be written in a similar way.The inductances LSfem and LRfem′, exactly their sum, has been determined indirectly by calculating the corresponding stator leakage linkage flux 1SL under the locked rotor simulation condition.The stator linkage flux 1S is generally defined for 2D problems by the formula: where As is the stator slot cross-section, Asi is the stator slot area of the all stator coils connected to one phase, A1 and A2 are the magnetic vector potentials at the beginning and the end of the turn side of the stator coils.Then, for the sum of the leakage inductances holds: where i1SL is the stator current under the locked rotor  The stator linkage flux Ψ 1S is generally defined for 2D problems by the formula: where A s is the stator slot cross-section, A si is the stator slot area of the all stator coils connected to one phase, A 1 and A 2 are the magnetic vector potentials at the beginning and the end of the turn side of the stator c 2018 ADVANCES IN ELECTRICAL AND ELECTRONIC ENGINEERING coils.Then, for the sum of the leakage inductances holds: where i 1SL is the stator current under the locked rotor simulation condition.The calculated sum of the leakage inductances L σSf em and L σRf em has been split into two halves in the usual way.According to [21], the difference between the magnetic vector potential in the middle of the stator slot and the air gap is very small under the no-load simulation condition, so the magnetizing inductance can be estimated as: where Ψ 1So is the stator linkage flux and i 1So is the stator current under the no-load simulation condition.
The rotor resistance R 2 referred to the stator has been calculated from the rotor power loss, namely from the rotor ring-end power loss P Rend , the outside rotor bar power loss P Bout and the rotor bar power loss P Rf em : where i R is the rotor current referred to the stator.Equations Eq. ( 17) to Eq. ( 20) hold at each time and have been enumerated over the period of 20 ms in the quasi-steady state using time-stepping analysis [22], i.e. for RMS values.For comparison, the calculated values of the motor parameters have been shown in Tab. 5 together with the ones calculated from the motor data plate and from the measured values.

Conclusion
The simulation results of the induction motor model match the motor catalogue values and the measured values rather well, especially for the nominal operating condition.There are only minor differences between the simulated and measured values of the stator current and motor torque, as seen in Tab. 2 for the nominal operating condition.In Tab. 3, the simulated currents of the rotor bar and rotor ring are compared to the currents from the theoretical presumptions.As seen, they match rather well.As regards the values of the stator current from the simulation under locked rotor condition, the difference is higher in comparison to the catalogue value, but it should be mentioned that the tolerance for the locked rotor current is stated as +20 % in the motor catalogue.Similarly, there is a certain difference between simulated and measured value of the stator current under no-load condition made by manufacturer as mentioned in Sec. 4. As seen in Tab. 5, the values of the equivalent circuit parameters of the tested motor resulting from its simulation are very close to their values obtained from measured parameters and from data plate parameters, which are listed in Tab. 1.
It is worth noting that the values of the simulated parameters depend on the accuracy of identification and calculation of the induction motor parameters used in the model, on the machine geometry (e.g.manufacture's stator and rotor laminations are necessary) and material properties, on the proper meshing especially of the motor air gap, and on the adjusted time step of simulation.
In the future, this model will be used for some studies of the tested motor under electrical and mechanical faults or anomalies.

Fig. 4 :Fig. 5 :
Fig. 4: Waveform of the current in the rotor ring at motor nominal speed.

Fig. 11 :
Fig. 11: Magnetic flux density and equipotential lines of magnetic vector potential under locked rotor at the nominal supply voltage.

Tab. 4 :
Comparison of the stator current under locked rotor.

Fig. 16 :
Fig. 16: Equivalent phase circuit of the induction motor. ) Fig. 7: Waveform of the inner torque and its mean value at the nominal speed.Tab.2: Comparison of the stator current and motor torque for the nominal operating condition.
Tab. 4: Comparison of the stator current under locked rotor.
Tab. 5: Comparison of the induction motor parameters.