An extended winding function model for induction machine modelling considering saturation effect

Winding function (WF) method is one of the electric machines modelling methods, which is used to simulate electric machines under different conditions. It has less simulation time and memory than other methods and, at the same time, has comparable accuracy over other methods. One of the significant challenges of WF ‐ based models is how to model saturation. The accuracy of the fault detection method in electric machines depends on the accuracy of the model and the method used for the simulation. A new extended model of the WF is presented for saturation modelling, which considers all parts of the core including teeth, slots, and yokes. In this method, the effect of the iron parts is emulated as a partial value that is added to the length of the actual air gap. The simulation results for normal and faulty conditions are shown in two cases, with and without considering saturation effect. Then the importance of saturation modelling has been demonstrated. Simulation is done in MATLAB software. Moreover, the experimental results are shown to verify the proposed model.


| INTRODUCTION
Modelling of electric machines in order to designing, fault detection, and control purposes have always been done in different ways. The two main features considered to compare different methods are the accuracy and complexity of the simulation model and runtime. Obviously, in some cases, the precision and complexity of the approach are essential, while in others, simulation time is more urgent. For example, for designing purposes where simulation time is not vital, more precise models can be used; however, for control purposes or identifying online phenomena, simulation time is important. Therefore, techniques that require less simulation time and maintain the simulation accuracy as much as possible can be suitable for many applications [1,2].
In the analysis of electrical machines, analytical methods using the mathematical transform models with many assumptions are unable to accurately analyse the electrical machines, especially under fault conditions [3][4][5]. In these analyses, the machine's air gap length is assumed to be uniform, and the spatial harmonics derived from the winding distribution and non-uniform air gap are neglected.
In contrast, magnetic field analysis approaches with acceptable accuracy have very low speed and very high memory requirements [6,7]. Equivalent magnetic circuits schemes, another method used for modelling induction machines, is based on the definition of equivalent magnetic reluctances. It defines a reluctance for each of the element of the machine; thus, defining many magnetic and electric circuits results in model complexity and high computational time [8,9]. Another method uses steady-state equivalent circuits such as T, Π, Γ, etc., which are simple, but these circuits are not suitable for transient state analyses [10].
Another analytical method that simulates a machine in the real space and is used in simulating electric machines under different conditions is called the winding function (WF) approach. This method analyses the machine based on the calculation of the machine inductances and their use in forming the electromagnetic coupling equations between the rotor and the stator [11][12][13][14][15][16]. Since there is no constraint on the stator winding and rotor bar positioning in the WF method, all the spatial harmonics obtained from the winding distribution and slots are included in the modelling. This approach can model many phenomena and has less simulation time and memory than other methods. At the same time, it has comparable accuracy to other methods, and can easily be used to analyse the transient state of electric machines under various conditions. Toliyat and Lipo, obtained the shape of the inductance functions of the squirrel cage machine by taking into account the winding distribution and the WF scheme. Then they wrote the differential equations governing the machine and simulated the transient and steady behaviour of the machine [11].
Moreover, winding faults modelling and air gap nonuniformity faults analysis due to static rotor eccentricity in induction machines have been performed using the WF method [17][18][19]. The effect of slot skewing and linear increase in magneto motive force (mmf) across the slots are modelled, and the impact of these factors on self and mutual inductances are investigated [20]. Besides, the WF approach is used to simulate the dynamic eccentricity of the rotor in the induction machines [21][22][23]. A developed model of WF was proposed for fault analysis of the induction machine under saturation condition. It modelled the saturation effect by modifying the length of air-gap and introduced a new way to find the air-gap flux density position [13]. In 2008, Ojaghi expanded on the previous work to incorporate the saturation effect into modelling, and proposed another way of finding the flux density of the air gap and modelling the saturated induction machine [14]. The tooth effect is defined as a coefficient in terms of the ratio of saturated and unsaturated voltages, which is added to the air-gap length. The calculation of this coefficient requires voltage measuring and using a lookup table.
In most of the simulations done so far using the WF, either the effect of the iron part has been ignored, or only the tooth and slot effects have been taken into account. The magnetic flux density, especially in transient states, is very high in the tooth areas, and in these areas, the core is saturated. If saturation is not considered in these areas, the resultant flux density is much higher than the saturation flux density. As a result, the resultant inductances are much higher than the actual value, which in turn increase the simulation error.
Considering the saturation effect for modelling and diagnosing fault types is much more critical. The efficiency of the fault detection approach in electric machines depends on the accuracy of the model and the method used for the simulation. From the technical and economic justification, the most suitable method for fault detection is the stator current analysis method. With the assumption of linearity of the core, the frequency components induced by the saturation in the stator current spectrum are not seen in the simulation signal, and the simulation signal spectrum differs from the real signal. Therefore, the implementation of fault detection systems without considering the saturation effect in modelling has errors.
A new extended model of the WF is presented for saturation modelling, which considers all parts of the core including tooth, slots, and yokes. In this method, the effect of the iron parts is added as a corrective value to the length of the actual air gap, which is obtained by an iterative algorithm. Moreover, for each point of the circumference, a corrective value is obtained.
In the remainder of the article, we will briefly refer to the WF method in Section 2. Then, in Section 3, we present the proposed model of WF. The fourth section presents the simulation results. The experimental results are shown in Section 5. Finally, the article concludes.

| WINDING FUNCTION APPROACH
In WF approach, the inductance relation of the windings is defined as following. In this approach, all factors affecting the non-uniform air gap, including the presence of slot and tooth, the asymmetry of the air gap due to rotor bar skewing, rotor eccentricity are considered [24].
F I G U R E 1 Turn function and winding function for a four-pole winding where, N A is the WF, n B is the turn function, r is the mean radius, and g is the air-gap length, which all depend on the angular position of the rotor, the angular reference on the stator, and the longitudinal position on the machine axis. Therefore, this method can calculate induction machine inductances under various asymmetries. The WF and turn function are shown in Figure 1 for a four poles winding.
Calculation of the air gap length in the presence of slot and tooth, taking into account the fringing effect, is usually based on the Carter model [25] illustrated in Figure 2.
In this case, the air gap is divided into three parts, stator side, middle part, and rotor side.
The effective air-gap length of the stator/rotor side is obtained as [25]: In which, g equals to one third of the air-gap length, g ¼ g 0 =3, β, and o 0 are defined as below: ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi In which, o is the slot opening. The total air-gap length considering the stator and the rotor slots is defined as below: In this model, the effect of the iron parts and, consequently, the effect of saturation have not been investigated. The next section presents a new modified model to take into account the saturation effect on the ferrous parts of the core.

| PROPOSED MODEL
Resultant mmf in the air-gap at any position, of the peripheral angle of the stator (φ) and the rotor angle (θ r ), is the sum of the mmf resulting from the stator circuits and the rotor circuits: The resultant mmf is dropped across the air gap and the iron parts: In previous studies, in order to introduce saturation effect, air-gap conductivity has been proposed as a function of the saturation coefficient, which is the ratio of the main components of the air-gap voltage under saturation and unsaturation conditions [10,11]. And the air-gap is assumed to be a function of the position and the flux density level of the air-gap, which these coefficients are approximate values, so the accuracy of these models is not high.
We introduce the iron parts directly into the problem and finally add their effect as a partial correction to the air-gap length of unsaturated condition.
In the proposed model, the mmf drop across the iron parts is modelled by changing the air-gap as below: where, Δgðθ; φ r Þ is the partial air-gap and g new ðφ; θ r Þ is the modified air-gap length. Calculation of the equivalent partial air-gap Δgðθ; φ r Þ, which is a different value for each operation point, is obtained by an iterative algorithm, whose implementation process is as follows.
First, we analyse the induction machine by assuming that it is linear (ignoring the drop of mmf across the iron parts). Then we obtain the total mmf using Equation (3), and the mmf drops across the air-gap and the iron parts of the core are obtained using Equation (4). However, to analyse the relationship Equation (4), a process must be performed as follows.
To solve Equation (4) and obtain mmf drop across each segment, it is necessary to define the mmf drops in terms of the magnetic field strength. The magnetic field strength in radial direction in the air-gap is considered constant, for simplicity. Therefore, the relationship of the mmf drop across the air-gap will remain as Equation (5).
On the other hand, the mmf drop across the iron parts is equal to the drop across the tooth plus the drop across the rotor and stator yoke.
where, the tooth drop is equal to the sum of the drops in the rotor and stator teeth.
The exact concept of the length l ths ; l thr is shown in Figure 3.
According to Figure 3, the following equations can be obtained for l ths ; l thr .
where, g ef f r ðφ; θ r Þ and g ef f s ðφ; θ r Þ are the effective air-gap functions obtained from the Carter relationship Equation (2) for the rotor and the stator. d ss and d sr are the depth of stator and rotor slots, respectively. The Carter relationship presents the effective air-gap function of a slotted surface versus a flat surface. Since in the induction machine, both the stator and rotor surfaces are slotted, therefore, we have to choose a hypothetical flat surface so that the air gap is divided into three parts, and the middle part acts as a flat surface for the two stator and rotor surfaces. Also, the yoke drop is equal to the sum of the drops in the rotor and stator yokes.
And the drop in stator and rotor yokes is obtained from the following equations.  where, r ms and r mr are the mean radius of the rotor and stator yokes. So far, the relationships are expressed based on the magnetic field intensity. To solve the equations, we need to express the relations based on the magnetic flux density.
The relationship between the magnetic flux density and the magnetic field strength over the air-gap and the iron parts is as follows.
where, G is a non-linear function obtained from magnetization curve of iron.
We now express the magnetic flux density in different parts in terms of the magnetic flux density of the air-gap.
The relationship of the magnetic flux density in the rotor and stator teeth will be as follows.
In these equations, the air-gap is considered to be the effective air-gap, so we consider the flux density of the air-gap to be equal to the flux density of the teeth.
And the magnetic flux density in the yokes is obtained as follows.
The A ys and A yr , respectively, are the cross-section of the yokes of the rotor and the stator in the circumferential direction. The flux at different points of the rotor and stator yoke is obtained from the following relationships.
The r s and r r , respectively, are the inner radius of the stator and the outer radius of the rotor, and l is the effective length of the machine.

F I G U R E 4 Flux density direction over the air-gap and yokes GHOLAMI ET AL.
φ max is an angle where the flux density of the air gap is maximum in absolute terms, noting that as many poles as a machine have, the maximum points we have. Usually, the maximum points are in the middle of the poles. The flux density vector inside the yoke is in the circumferential direction, depending on the location of the yoke periphery, whether it is with the N or S poles, can be triangular or countertriangular. This principle is illustrated in Figure 4.
Thus, the right-hand components of the relationship Equation (4) are written in terms of the magnetic flux density. It can be seen that this relation is a non-linear equation so that numerical methods can be used to solve it. The Newton-Raphson iteration method is used. Therefore, the mmf drops across the air-gap, and the iron parts can be calculated. Now the equivalent air-gap is calculated using the Equations (6) and (7). Returning to the beginning of the simulation, we compute the inductance functions using Equation (1) with the new airgap and repeat the process until the difference of one of the quantities, such as machine currents, is less than the acceptable tolerance in two consecutive iterations.

| SIMULATION RESULTS
Simulations have been performed for two cases, healthy motor and four broken-bar motor for a 3 hp, 380-V, four-pole, threephase, star-connected squirrel cage induction motor. The simulation results are compared in two cases, with and without the saturation effect. Figure 5 shows the motor's starting currents at no load obtained from the simulation results in two conditions with and without the saturation effect. As can be seen, the F I G U R E 6 No-load starting speed; with and without considering saturation effect F I G U R E 5 Starting current in no-load condition; with and without saturation effect maximum starting current reaches 30.2 A with saturation. However, in the linear model, the maximum starting current reaches 23.5 A. Figure 6 shows the no-load starting speed variations in two modes with and without the saturation effect. As can be seen, given the saturation effect, the motor reaches its maximum speed in 35 ms. While ideally assuming the core, the motor reaches its maximum speed in 42 ms. Therefore, given the saturation, both the current and the starting torque increase.

| Healthy motor
Due to the high starting currents, the machine's mmf increases substantially, causing the core to saturate. When the motor is saturated, the magnetization inductance of the core decreases. Therefore, the magnetizing current and consequently, the motor current increases. As the current increases, the starting torque also does.
After the motor has reached steady state, a torque of 14 Nm is applied to it. Figure 7 shows the speed of the motor when it reaches steady state in two modes with and without the saturation effect. As can be seen, due to the saturation effect, the output speed decreases from 1429 to 1426 rpm, indicating a decrease in the saturation output power. Speed fluctuations also increase. Figure 8 also shows the motor current in this case. As can be seen from the saturation effect, the maximum value of the steady state current has increased from 7.1 to 7.8 A.

F I G U R E 8
Steady-state current at loading; with and without the saturation effect F I G U R E 7 Steady-state speed at loading; with and without the saturation effect Also, the frequency spectrum of the stator current in the main frequency band for two states with and without the saturation effect is shown in Figure 9. As observed, in the saturation state, in the main frequency band, the amplitude of the harmonics increases.
Another frequency band that is of interest in machine studies is the slot harmonic band [20], so the simulation results for the stator current spectrum in this band for two modes with and without saturation effect are shown in Figure 10. As can be seen, the principle slot harmonic (PSH) with no saturation effect (618) is 3 Hz higher than that of with saturation effect (615). The reason is the speed reduction considering the saturation effect. As shown in Ref. [20], the frequency of the main component of the slot decreases with decreasing speed and increasing slip. Also, in this band, components 245, 345, 515, 815, and 915 Hz are observed with significant values in the saturation mode, which are very small (or absent) in the no saturation mode.
Furthermore, in order to investigate the slot effect, the simulation results for stator current spectra in loading condition is shown in Figure 11. As can be seen, the slot has not significant effect on the main frequency band. However, the PSH band is almost eliminated.  Figure 12 shows the stator current spectrum in two modes without and with saturation effect, where the four rotor bars are broken. As can be seen in this figure, the breaking of the rotor bars has increased the amplitude of the lateral harmonics of the stator current. However, the increase in the harmonic components of the stator current in the case of without the saturation effect is more than the saturation effect case. The cause of this phenomenon can be stated in this way.

| Broken-bar fault
When the rotor bar is broken, the absence of current and opposite flux in the rotor bar increases the flux density around it. Therefore, breaking the rotor bar causes asymmetry of magnetic flux density in the air gap and, consequently, induction of frequency components in the stator current. If the ideal magnetic core is assumed, there is no limit to the increase in the flux density. Whereas, in nonideal core, an increase in the flux density in the broken bar region would result in local saturation of the rotor core. Therefore, given the saturation effect, the flux density is limited to the saturation density. In other words, the asymmetry of the flux density decreases. As a result, saturation reduces the fault-induced asymmetry and consequently reduces the amplitude of the harmonic components caused by the fault.

| EXPERIMENTAL RESULTS
A 3 hp, 380-V, four-pole, three-phase, star-connected squirrel cage induction motor, similar to the simulated one, shown in Figure 13, is used to experimental test. The experiment is performed for healthy conditions.
The starting current in no load condition is shown in Figure 14a. As can be seen the maximum current is around 30 A which confirms the simulation result in saturation mode in Figure 5.
A synchronous generator is coupled to the motor to perform the loading condition test. The loading current is shown in Figure 14b. The maximum value of the current is 7.9 A. This result confirms the simulation result in Figure 8.
The spectrums of the stator current in two bands, the rated frequency and the principle slot harmonic, are shown in F I G U R E 1 1 Stator current spectrum in case of without considering slot effect; (a) the main frequency band, (b) the PSH band. PSH, principle slot harmonic Figure 15a,b, respectively. The harmonics in the main frequency band have significant amplitude which confirm the saturation mode results in the proposed model ( Figure 9b). Moreover, the presence of harmonics in the PSH band is seen with significant amplitude which confirms the result of simulation in Figure 10b.

| CONCLUSION
Having a relatively accurate model of electrical machines requires consideration of various parameters affecting the performance of the machine. One of these effects, which has been less discussed in the WF approach, is the saturation effect of the core.
We propose an algorithm to consider the saturation effect in modelling and simulation of an induction machine by using the WF approach. In this algorithm, the machine is first analysed by assuming linearity, and the mmf drop in the air-gap and the iron parts is calculated. The mmf drop in the iron parts is replaced by a partial air-gap, and a new air gap function is obtained. Then by returning to the beginning of the simulation, the WFs and inductances are calculated with this new air-gap function, and the machine equations are reanalysed. This process is repeated until a parameter change reaches to an acceptable level. The simulation results for healthy and faulty motor are evaluated in two cases; with and without considering the saturation effect. The importance of considering the saturation effect on the results of the stator current spectra, used in fault detection systems, is shown. It changes the amplitude of the harmonics in the current spectra. Moreover, the experimental results are compared with the simulation ones. The simulation results are very close to the practical results, which confirm the accuracy of the proposed model. F I G U R E 1 2 Stator current spectrum for a four broken bars motor; (a) without the saturation, (b) with the saturation