Partial Inductance Model of Induction Machines for Fault Diagnosis

The development of advanced fault diagnostic systems for induction machines through the stator current requires accurate and fast models that can simulate the machine under faulty conditions, both in steady-state and in transient regime. These models are far more complex than the models used for healthy machines, because one of the effect of the faults is to change the winding configurations (broken bar faults, rotor asymmetries, and inter-turn short circuits) or the magnetic circuit (eccentricity and bearing faults). This produces a change of the self and mutual phase inductances, which induces in the stator currents the characteristic fault harmonics used to detect and to quantify the fault. The development of a machine model that can reflect these changes is a challenging task, which is addressed in this work with a novel approach, based on the concept of partial inductances. Instead of developing the machine model based on the phases’ coils, it is developed using the partial inductance of a single conductor, obtained through the magnetic vector potential, and combining the partial inductances of all the conductors with a fast Fourier transform for obtaining the phases’ inductances. The proposed method is validated using a commercial induction motor with forced broken bars.


Introduction
Induction machines (IMs) are key components of modern industrial processes, and their unexpected failures can generate unscheduled stops and high economic losses. To avoid this risk, on-line condition monitoring of IMs has become an integral part of industrial maintenance systems [1,2]. A great research effort has been devoted in the last decades to the development of fault diagnostic systems that can detect multiple machine faults in an early stage, giving a reliable information about the machine condition, and avoiding false alarms [3,4]. Different signals have been proposed to build fault diagnosis systems, such as vibrations [5], stator currents [6], thermal images [7], acoustic signals [8,9], etc. In particular, diagnostic algorithms that can detect the faults through their characteristic signature in the stator currents, known as motor current signal analysis (MCSA), have gain a great relevance [6,10], because they can be applied on-line using a simple hardware and software, just current clamps to measure the stator currents and the fast Fourier transform (FFT) to analyze them in the frequency domain [11,12]. Multiple types of faults can be detected by identifying the characteristic harmonics that they generate in the stator current, with frequencies that have been established theoretically.
and for mechanical looseness f loose [23] f loose = f 1 1 + k 1 − s p · n k = 1, 2, . . . , n = 2, 3, . . . (4) where s is the slip, f 1 is the frequency of the power supply and p is the pole pairs number. Other types of faults, such as bearing faults [24][25][26], or gearbox faults [27] generate components in the stator current with frequency expressions similar to Equations (1)- (4). These equations are valid both in steady-state regime and in transient regime, but in this case the frequencies given by Equations (1)-(4) change along the time [15], and advanced time-frequency transforms must be used instead of the FFT [25,28]. Commonly, fault diagnostic systems that are based on the analysis of the current signal, such as MCSA, try to measure the amplitude of the harmonic components with frequencies given by Equations (1)- (4) in the spectrum of the stator current. These values are then compared with predefined threshold alarms, in data based diagnostic systems [5], or with the output of machine models, in model based diagnostic systems [5,20,[29][30][31]. In both cases, accurate machine models are needed, for fixing the values of the threshold alarms, for running the model-based fault diagnostic system [32], or for developing new diagnostic algorithms [33]. The most accurate machine models are built using the finite elements analysis (FEA) [34], because it can take into account the actual geometry of the machine, including the asymmetries induced by the fault (eccentric machines), the characteristic of the magnetic materials, which depend on their operating point (saturation, hysteresis), and the configuration of the windings, which can be highly asymmetrical in case of rotor asymmetries or inter-turn short-circuits faults [35]. However, FEA models demand a great amount of computational resources and computing time [17]. This is a serious drawback for using FEA models on embedded on-line diagnostic hardware [24], such as DSPs or FPGAs [36], which have limited computing resources, or in systems which require a fast response, such as real-time model-based diagnostic systems [37]. A growing trend is the use of the IM model in hardware-in-the-loop (HIL) systems, which are used for developing and testing real-time diagnostic hardware, [38], or for training advanced artificial intelligence (AI) based diagnostic systems, such as the set membership identification (SMI) approach [14], neural networks [39][40][41], support vector machines [42] or fuzzy inference systems [43].
To alleviate these drawbacks, several solutions have been proposed in the technical literature, such the use of coupled finite elements and circuit equations [44], improved conformal mappings [45], or field reconstruction methods (FRM), based on linear superposition of FEA solutions [46]. An alternative approach to FEA models are fast and simpler analytical models of the induction machine [47], such as the magnetic equivalent circuit (MEC) model [36,48], which represents the induction machine as a mesh of reluctance elements, or the permeance network model (PNM) [49]. Recent works also suggest the use of generalized two-phase model of an induction machine of fifth order [50], space phasor models [27], complex-vector models [51], or equivalent circuit models [52].
Among the analytical models, the multiple coupled circuit model (MCCM) [53,54], also called the winding function approach [55], has been extensively used for modeling the IM in faulty conditions [56][57][58]. As Ojaghi and Faiz [55] pointed out, it is the most detailed and complete model used to analyze the performance of cage IMs, especially under various faults. It relies on the time-varying self and mutual inductances of the machine phases and rotor-loops, and their derivatives, which are calculated by integral expressions representing the placement of the winding turns along the air-gap periphery. However, these integral expressions must be computed for every rotor position and phase winding layout, which can be altered due to the asymmetries generated by the machine's fault, resulting in the need of heavy computing resources to build the machine inductances matrices under faulty conditions. As it is stated in [54], this task typically consumes a high amount of time, so that only discrete curves of inductance versus rotor position are calculated and linear interpolation is applied at intermediate rotor positions. To solve this problem, in this paper, a novel, fast approach is proposed for obtaining the matrices of self and mutual inductances, making use of the concept of partial inductance of a single conductor [59,60]. Instead of solving the integral expression for every different winding layout, the partial inductance of a single conductor [61] is computed, using the magnetic vector potential (MVP), and the total self and mutual inductances of the windings are obtained by simply combining the partial inductances of their conductors. This combination is made in a very fast way using the fast Fourier transform. The partial inductance approach has been used in the partial element equivalent circuit (PEEC) method [62] for the electromagnetic modeling of power electronic modules [63], antennas [64], high voltage transmission lines [65] or printed circuit boards (PCBs) [66], but it is the first time, to the authors' best knowledge, that it is proposed for obtaining the phase inductances of IMs. Compared with the WFA, this novel method is characterized by the following main points: 1.
The conductor, instead of the coil, is used as the basic winding unit, which simplifies the modeling of arbitrarily complex winding layouts.

2.
The partial inductance of a single conductor, computed using the MVP, is taken as the characteristic function of the machine's windings, instead of the winding functions based on the magnetomotive force (MMF) of the coils.

3.
A fast circular convolution, computed through the FFT, is used to obtain the mutual and self inductances of the phases, instead of solving the integral expressions of the winding functions for every rotor position.
The structure of this paper is as follows. In Section 2, the electro-mechanical equations used for modeling the IM are briefly presented. In Section 3, the partial inductance of a single conductor is obtained through its MVP, and in Section 4 the FFT is used for assembling the matrices of self and mutual inductances of the machine, combining the partial inductances of the phases's conductors. Section 5 presents the experimental validation of the proposed approach, using a commercial IM with broken bar faults. Finally, Section 6 presents the conclusions of this work.

Electro-Mechanical System of Equations of an Induction Machine
The following system of equations can be written for an induction machine with m stator and n rotor phases with arbitrary layout, that is, even under fault conditions, [ where T e is the electromechanical torque of the machine, T L is the load torque, J is the total system inertia (rotor plus load), Ω is the mechanical speed, θ is the angular position of the rotor, and ϕ is the angular coordinate. Two additional constraints, namely the sum of all the stator currents equals 0, and the sum of all the rotor currents equals 0, must be added to this set of equations. To solve the system of Equations (5)- (12), the inductance matrices and their derivatives must be calculated very accurately, in order to be able to reproduce the fault harmonics given by Equations (1)-(4). The inductances in Equations (7), (8) and (11) correspond to the self and mutual inductances of the machine's phases, and of the rotor loops in the case of squirrel cage machines (two consecutive bars and the sections of the end-rings between them). The inductances of these windings, such as the ones represented in Figure 1, can be computed summing up the coil inductances for all the phases' coils, as in WFA. However, this approach must deal with the wide variety of coil types that can be used for building a phase winding. Besides, these inductances may change with the rotor position, as in the case of rotor-stator mutual inductances, and also in the case of winding faults, such as inter-turn short-circuits, or bar breakages. This approach results in complex and lengthy computations that typically consume a high amount of time [54]. The use of the partial inductance concept can greatly simplify this process, because the partial inductance of the conductor, instead of the coil inductance, is chosen as the building block for the computation of the inductance matrices, which is performed using a very simple and efficient FFT-based procedure, even in the case of complex or faulty winding layouts. In this work, the following assumptions will be made for calculating the phase inductances, common to many analytical IM models presented in the technical literature:
The eddy currents are neglected. The windings are considered to be made of filamentary conductors placed on the external rotor surface or on the internal stator surface. 4.
The end-coil leakage inductances need to be pre-calculated, and are treated as constants in Equations (7), (8) and (11), using explicit expressions that can be found in [52,68].
Some of these assumptions can be removed by using a modified inverse air-gap function which takes into account the saturation [55], the rotor eccentricity [30], or the slotting effects [69]. Nevertheless, as the focus of this paper is the introduction of the partial inductance model of the IM, a constant air-gap has been used to present the proposed approach, applied to the simulation and experimental test of an IM with a rotor asymmetry (broken bars fault), the same fault analyzed in [11,12,14,15,56], among many others.

Phase Inductance as a Sum of the Conductors' Partial Inductances
Let us consider a single coil of the machine windings, represented as a rectangular current loop in Figure 2, where conductors c 1 and c 2 are parallel to the machine's axis, and conductors c 3 and c 4 belong to the end windings. The self inductance of this single coil, under the assumptions given in Section 2, is given by Equation where Ψ and I are the flux linked by the coil and its current, respectively, assuming that I is the only current in the machine. The same equation holds for the mutual inductance between two coils, just using the current of the second coil in Equation (13). However, the coil inductance given by Equation (13) is a global coil characteristic which can't be assigned to any segment of the loop represented in Figure 2. On the contrary, a partial inductance can be assigned individually to each of the conductors that forms the coil loop. In effect, Equation (13) can be expressed as Equation where B is the magnetic flux density through the loop area S. Using B = ∇ × A, where A is the magnetic vector potential, and applying Stokes' theorem, gives Equation where C is the contour of the current loop. However, the contour integral in Equation (15) can be broken into the four segments of the loop of Figure 2 as where C i is the contour of each segment of the loop. The net inductance of each segment C i , L i in Equation (16), represents the contribution of the segment C i to the total inductance of the current loop.
Using the concept of partial inductance, the inductance of each segment i of the coil in Equation (16) L coil i (i = 1 . . . 4) can be expressed as the sum of the partial inductances of segment i in Equation Lp ij (17) where the partial inductance Lp ij is defined [59,60] as the ratio between the magnetic flux crossing the rectangular surface between the segment i and infinity, and the current I j that produces that flux (which, in this case, is the same for all the loop conductors). That is, Using Equations (16) and (17), the total inductance of the rectangular loop in Equation (13) can be expressed as the sum of all the partial inductances of its segments, as As stated in the previous section, the end-coil leakage inductance is considered in this work as a constant value, computed a priori, so Equation (19) must be applied only to the axial conductors of the coil, C 1 and C 2 in Figure 2. With this assumption, the sums in Equation (19) must be extended only to these two conductors, as In a similar way, the self-inductance of a phase A with N A conductors can be computed extending Equation (20) to account for all the phase's conductors, as where L end_phase represents the total end-winding leakage inductance of the phase. The expression for the mutual inductance between phase A and other phase B, with N B conductors, is analogous to Equation (21), Lp ij (22) with the additional assumption of neglecting the mutual inductance between the end-coils of both phases.
Comparing Equations (21) and (22), it is clear that L phase A in Equation (21) can be obtained using Equation (22) as so this work will focus in the effective computation of Equation (22) The practical implementation of the proposed method seems to be cumbersome. First, a partial differential equation (PDE) must be solved for every conductor to obtain the MVP in Equation (18). Second, two double summations must be made in Equations (21) and (22), for all the stator and rotor phases. This process must be repeated for every position of the rotor, for obtaining the mutual inductance between a rotor phase and a stator phase. Nevertheless, in the next sections, an analytical expression for Equation (18) is derived, and a simple procedure is given for obtaining Equation (22), for any possible displacement between phases A and B, using a single FFT.

Partial Inductance between Axial Conductors
As only the partial inductances between the axial conductors of the phases are considered in Equations (21) and (22), only the axial component of the MVP must be taken into account. Neglecting end effects, a central cross section, perpendicular to the machine axis, is analyzed, thus using a 2D approximation for computing the MVP. Under these conditions, the MVP has a single, axial component, A z . If the length of the axial conductors is equal to the effective length of the magnetic core of the machine l maq , the value of Lp ij in Equation (18) is simply In an induction machine, the boundary condition at the infinity in Equation (18) can be replaced with a zero Dirichlet boundary condition at the external surface of the stator, where the MVP is assumed to be zero (no flux crossing the external surface of the stator). With this boundary condition, Az ij in Equation (18) is just the MVP at the conductor's position, which is equal to the yoke flux per unit length at this position.
The method used in this work for obtaining the partial inductance Lp ij in Equation (24) between two axial conductors i and j consists in feeding conductor j with a unit current, being this the only current in the machine, and computing the MVP at the position of conductor i, so that Equation (24) becomes For applying the proposed method, it is necessary to obtain the MVP generated by a single conductor, fed with a unit current. The MVP generated by a linear current in an air-gap with concentric circular cylindrical boundaries, as shown in Figure 3, must satisfy the Poisson's equation Equation (26) can be solved using the method of separation of variables, as in [70,71]. Placing the conductor at the origin of the angular coordinate (ϕ = 0 in Figure 3), and assuming that the iron of the cylindrical boundaries is infinitely permeable, the solution of Equation (26) for the conductor shown in Figure 3 is The constant, cyclic part of Equation (27) (ln(c), ln(r)) can be discarded because, in any induction machine, the net sum of the linear currents through a transversal section is always zero. Figure 4 plots the harmonic part of Equation (27) for an example machine with a = 1, c = 1.4 and b = 1.5. If the conductors are assumed to lie on the outer rotor surface or in the inner stator surface ( Figure 5), then the linear currents are restricted to c = a or c = b in Figure 5. Besides, for obtaining the partial inductances of the conductors, the MVP must be calculated only on these same surfaces, that is, r = a or r = b in Figure 5. These four possibilities (c = a, c = b, r = a, r = b) reduce Equation (27) to only two different expressions: the MVP generated by an axial conductor on the same surface where it lies, or on the opposite one: The partial inductance of Equation (25) between two axial conductors i and j can now be computed using the values given by Equation (28). If conductor i is placed at the origin (ϕ = 0 in Figure 3), and conductor j is placed at an arbitrary angular position ϕ, then the partial inductance between these two axial conductors is If now conductor i is assumed to be at an angular position ϕ = α, instead of ϕ = 0, then the partial inductance with conductor j will be simply The derivative of the partial inductance with respect to the angle, which is used for computing the torque in Equation (11), can be obtained directly from Equation (29) as In Equation (31), only the derivative of the mutual inductance between conductors in opposite surfaces is needed, due to the uniformity of the air-gap.

Discretization of the Expression of the Partial Inductance between Axial Conductors
Equation (29) must be discretized in order to be used in an induction machine's model. Therefore, the air-gap circumference is divided into N equal-length intervals, with an angular width equal to ∆ϕ = 2π/N. In this case, the maximum number of harmonics of the infinite series in Equation (29) that are needed to represent the partial inductance is, by Shannon's theorem, N/2. Thus, the sequence that gives the partial inductance between a conductor placed at the origin, and other conductor j placed at an angular position j · ∆ϕ, with j = 0 . . . N − 1, is Analogously to Equation (30), if conductor i is not placed at the origin, but at an angular position i · ∆ϕ, then Equation (32) becomes where ((j − i)) N = |j − i| mod N is understood as a modulo N operation. In this work, to simplify the mathematical expressions, all the matrix indexes are considered as modulo N, without using the full notation of Equation (33). Following this convention, Equation (33) is expressed as The derivative with respect to the angular coordinate of the partial inductance between two conductors placed on opposite surfaces, given by Equation (31), can be also discretized with the same procedure, giving a sequence defined by

Vectors Containing the Distribution of Conductors of Phases A and B and the Partial Inductances between Two Axial Conductors
The process for computing the mutual inductance between two phases A and B, such as those represented in Figure 5, is based on the substitution in Equation (22) of Equation (33). To perform this process, using discrete arithmetic, three column vectors of N elements are needed.

•
The vector L p0 , whose element j, L p0 [j], contains the partial inductance between a conductor placed at the origin and other conductor placed at an angular position j · ∆ϕ, using the corresponding Equation (32) (depending if both conductors are in the same surface or in opposite surfaces).
• The vector dL p0 , whose element j, dL p0 [j], contains the derivative with respect to the angular coordinate of the partial inductance between a conductor placed at the origin and other conductor placed at an angular position j · ∆ϕ, using the corresponding Equation (35).
• The vector Z A0 , whose element j, n A0 [j], contains the number of conductors of phase A located in the angular interval of length ∆ϕ centered at position j · ∆ϕ. The direction of the currents at each position is represented using positive and negative values, depending on the current direction In the definition of Z A0 , it is assumed that the axis of the distribution of conductors of the phase coincides with the origin of coordinates ϕ = 0. If the phase is shifted by an angle ϕ A = k A · ∆ϕ, then all the elements of Z A0 are shifted by k A , that is • The vector Z B0 , which contains the distribution in the airgap of the conductors of phase B, which is defined in the same way as vector Z A0 In a similar way to the definition of Z A0 , it is assumed that the axis of the distribution of conductors Z B0 coincides with the origin of coordinates ϕ = 0. If the phase is shifted by an angle ϕ B = k B · ∆ϕ, then all the elements of Z B0 are shifted by k B , that is It is worth mentioning that no restriction has been put on the definition of vectors Z A0 and Z B0 , so that arbitrarily complex winding layouts can be used in the proposed model, including the uniform distribution of the conductors along the slot width, or along skewed slots.

Vector L AB Containing the Mutual Inductances between Phases A and B for All Their Relative Positions
With the vectors defined in Section 4.1, the expression of the mutual inductance between phase A, with its axis located at an angular position ϕ A = k A · ∆ϕ, and phase B, placed at other angular position ϕ B = k B · ∆ϕ, L k A k B , can be formulated, using Equation (22), as Equation (42) can be formulated with vectors L p0 , Z A0 , Z B0 , using Equations (34), (39), and (41) respectively, as Using the change of variables j = j − k A , Equation (43) becomes and, using the change of variables i = i − k A , Equation (44) becomes That is, the mutual inductance between phases A and B depends only on their relative position (k B − k A ) · ∆ϕ. In this way, it can be defined a vector L AB whose element k, L AB [k], contains the mutual inductance between phases A and B for an angular separation between them k · ∆ϕ Equation (46) is simple to setup, because it contains only the expression of the partial inductance between two axial conductors (Equation (32)), and the distribution of the conductors (not the coils) of phases A and B, all of them obtained using a common origin of coordinates. However, at the same time, it is cumbersome to calculate, because a double summation is needed for obtaining the mutual inductance L AB [k], for each value of k at N different positions. Nevertheless, this process can be performed in an extremely efficient way using the FFT, as proposed in the present work.

Matrix Formulation of the Expression of the Mutual Inductance between Phases A and B
Equation (45) can be expressed in matrix form as where TcL p0 is a square N × N matrix defined as and TcZ B0 T is the transpose of a square N × N matrix defined as Matrix TcL p0 is built using the column vector L p0 in Equation (36) as its first column, and rotating each successive column by one position. That is, TcL p0 is a Toeplitz circulant matrix. Matrix TcZ B0 has been built using the same procedure, thus it is also a Toeplitz circulant matrix.

Computation of the Mutual Inductance between Phases A and B Using the FFT
A nice property of the Toeplitz circulant matrices is that they become diagonal by a transformation into the frequency domain. Therefore, Equation (48) can be reduced to a simple element-wise multiplication of three vectors in the frequency domain, using the FFT (and its inverse, the IFFT), as where the superscript * stands for the conjugate of each vector element, and the operator • stands for the element-wise or Hadamard product of two vectors (written as "*" in MatLab environment). This equation is similar to the one presented in [72], where the convolution theorem has been applied to obtain the mutual inductances between two phases, based on the magnetomotive force generated by a single conductor, instead of the MVP. Equation (50) gives, with a single element-wise multiplication of three vectors in the frequency domain, the values of the mutual inductance between phases A and B for every angular displacement between them (with a precision ∆ϕ). It relies only in the distribution of the phases' conductors, and on the analytical expression of the partial inductance between two axial conductors, Equation (32), which makes it very simple to implement. This approach differs from the approaches followed in [70,71], where the equation of the MVP in cylindrical coordinates (Equation (26)) must be solved again for each distribution of phase conductors, which limits its practical application to simple distributions, such as the sinusoidal ones.
The derivative with respect to the angular coordinate of the mutual inductance between phases A and B, located on opposite surfaces, has the same expression as Equation (50), just replacing the sequence L p0 in Equation (32) by its angular derivative, the sequence dL p0 in Equation (35), giving The procedure for obtaining the mutual inductance between two phases A and B, L AB in Equation (50), has been depicted graphically in the block diagram shown in Figure 6.
The same block diagram in Figure 6 is valid for computing the derivative of the mutual inductance between phases A and B, dL AB in Equation (51), just replacing in Figure 6 the vector L p0 in Equation (36) by its derivative, the vector dL p0 in Equation (37).
Amgular separation between phase A and phase B

Experimental Validation
The proposed method has been applied to the fault diagnosis of a commercial squirrel-cage induction motor, whose characteristics are given in Appendix A. The type of fault that has been used for the experimental validation of the proposed model is a broken bars fault. The use of the proposed partial inductance approach is particularly well suited for modeling a squirrel cage induction motor with broken bars. In this case, instead of using rotor meshes formed by two consecutive bars, as usually done, the faulty rotor winding can be modeled using directly the bar currents in Equation (10), and the partial inductances of the rotor bars in Equations (7), (8) and (11). Following this approach, the presence of multiple broken bars, consecutive or not, can be introduced in the model just by eliminating the broken bars entries in the current and voltage vectors, Equations (9) and (10) , instead of modifying the inductances of the rotor meshes affected by the broken bars. This approach avoids the cumbersome process of recomputing the rotor inductances matrix to take into account the variations in some rotor meshes generated by the broken bars fault.
It is worth mentioning that the election of a commercial squirrel cage motor for performing the experimental validation of the proposed method does not exclude its validity for other types of induction machines, such as wound rotor induction machines. In fact, Equations (50) and (51) do not impose any constraint on the winding distribution, and are equally valid for both single rotor bars and wounded rotor phases.

Experimental Setup
The test equipment, displayed in Figure 7, consists of a current clamp, whose characteristics are given in Appendix B, a 200 pulse/revolution incremental encoder, a Yokogawa DL750 Oscilloscope and a Personal Computer connected to it via an intranet network. As shown in Figure 7, the broken bar fault has been artificially produced by drilling a hole in the selected bars. The experimental validation of the proposed approach has been carried out using a set of artificially damaged rotors with one or two broken bars, beginning with two consecutive broken bars (Positions 1-2), and increasing gradually their separation (Positions 1-3, 1-4, etc.), up to a separation of seven slots between the two broken bars (Positions 1-8), as shown in Figure 8. An additional healthy rotor, with no broken bars, has also been used for comparison purposes. This gives a total number of nine experimental tests: a healthy motor, a motor with one broken bar, and seven motors with two broken bars at different positions.
The same stator has been used in all the experimental tests (Figure 9), which allows comparing the effects of the broken bar fault in the stator currents in a controlled environment. To this end, the motors were disassembled, and one of the stators was mounted in the test bed. All tests have been carried out mounting the various rotors (see Figure 9) in this stator. The induction motor under test (Appendix A) is connected via a belt coupling to a DC generator, which feds a resistive load (see Figure 10). Both the resistive load and the field excitation of the generator can be controlled, so that the induction machine works at rated speed 1410 r/min (s = 0.06). The current of a stator phase has been measured using a sampling frequency of 5000 Hz, during an acquisition time of 50 s, and its spectrum has been obtained for identifying the fault harmonics given by Equation (1).   The induction motor under test (Appendix A) is connected to a DC generator via a belt coupling. The DC machine feeds a resistive load. Both the resistive load and the field excitation can be controlled so that the induction machine works at rated speed.

Model Setup
The partial inductance between two conductors of the simulated machine in Equation (32) has been represented in Figure 11, for the cases of the two conductors lying in the same surface (Figure 11, top), or in opposite surfaces (Figure 11, middle). The derivative of the mutual inductance of a stator conductor and a rotor one in Equation (35) has been also represented in Figure 11, bottom. These values depend only on the geometrical dimensions of the machine. The distribution of the conductors of a stator phase and the one that corresponds to a skewed rotor bar of the simulated machine are shown in Figure 12. A direct substitution in Equations (50) and (51) of the values presented in Figures 11 and 12 gives the mutual inductance between two phases, for every relative angular position, and its derivative respect to the angular coordinate. The mutual inductance between a stator phase and a rotor bar of the simulated machine, as a function of the rotor position, is given in Figure 13, top, and its angular derivative is given in Figure 13, bottom. The simulations with the motor model have been carried using a constant load torque T L in Equation (12), with a value equal to the motor rated torque, that is, T L = 9550 × 1.1/1410 = 7.45 Nm.

Diagnosis of a Single Broken Bar Fault
The procedure for the diagnosis of a single broken bar fault is based on the MCSA technique. In the case of a broken bar fault, the magnitude of the characteristic fault harmonics given by Equation (1) increases significantly. For the tested and simulated motors, the rotor speed is the rated one (1410 r/min , f 1 = 50 Hz). The main fault harmonics used for the diagnosis are those corresponding to a value k = ±1 in Equation (1): the lower sideband harmonic (LSH), with a frequency f LSH = (1 − 2s) f 1 , and the upper sideband harmonic (USH), with a frequency f USH = (1 + 2s) f 1 . In the case of the tested and simulated motors, s = 0.06, so that f LSH = (1 − 2 × 0.06) × 50 = 44 Hz, and f USH = (1 + 2 × 0.06) × 50 = 56 Hz. Figure 14 shows the spectrum of the currents measured in the experimental tests, with a healthy and a faulty motor with a single broken bar. Figure 15 shows the same spectra obtained from the simulated motor.
The magnitude of the fault harmonics are presented in Table 1, showing a good agreement between the simulated and the experimental data. It is worth mentioning that, in the case of the experimental motor, two additional harmonics appear at frequencies of 43.5 Hz and 56.5 Hz, due to the belt used for coupling the load to the test bed. In fact, when the motor is tested unloaded, with the belt removed, these harmonic do not appear, so they are probably generated by an axial eccentricity induced by the asymmetric load coupling to the motor shaft. These harmonics do not appear in the case of the simulated motor, because the partial inductances model presented in this work does not take into account the effect of axial eccentricity. 40       Major motor manufacturers have reported cases where the damaged bars appear randomly distributed around the rotor perimeter, indicating that the failure of non-adjacent bars is fairly common in large cage induction motors. Figure 16 shows the rotor of a 3.15 MW, 6 kV induction motor with a double breakage fault, affecting to non-consecutive rotor bars. In the case of a double breakage, the magnitude of the LSH is a function of the relative position of the two broken bars. In [73] it has been shown that the ratio between the LSH in the case of double and single breakages depends on the angle between the broken bars as Equation where p is the number of pole pairs and α bb is the angle between the two broken bars. From Equation (52), it can be deduced that if α bb approximates π/2p, that is, half a pole pitch, then the second breakage reduces the magnitude of the LSH to a value lower than its magnitude in the case of a single breakage. Therefore, in this case, a motor with two broken bars could be erroneously diagnosed as a healthy motor. This behavior is more challenging to simulate than the single broken bar fault, and it has been selected to verify the validity of the proposed model, based on partial inductances, for fault diagnosis.
The following procedure has been followed to perform the experimental tests: 1. Seven rotors in healthy condition have been successively mounted and tested, to verify that there are no significant differences between the tested rotors.

2.
A bar has been drilled in each of the rotors, denoted as b 1 bar. The motors have been tested at rated speed, and, using the spectrum of one of the stator phase currents, the magnitude of the LSH (LSH single in Equation (52)) has been recorded.

3.
A second bar has been drilled at different positions from the first one (see Figure 8). The motors with a double breakage have been tested again at rated speed, and, using the spectrum of one of the stator phase currents, the magnitude of the LSH (LSH double in Equation (52)) has been recorded.

4.
Finally, the ratio LSH pu has been computed for every rotor, and compared with the same value obtained with the simulated motor. The results are presented in Figure 17 and in Table 2. The time spent in each simulation (50 s in steady state) was 12 s, using a computer whose characteristic are given in Appendix C.  Figure 17. Comparison between the simulated and the experimental data obtained for the motor with one broken bar (first column), and with two broken bars separated a variable number of healthy bars. The results shown in Figure 17 and in Table 2 indicate that the results obtained with the proposed model clearly follow the experimental trend, which confirms its validity as a tool for fault diagnosis. It is worth mentioning that, with the proposed approach, all the experimental tests, using the rotors presented in Figure 8, have been simulated using a single inductances matrix, which contains the partial inductances between stator phases and rotor bars. Each type of the considered faults has been simulated just by eliminating the faulty bars from the set of Equations (7), (8) and (11), without any modification to the rest of the terms in these equations.

Conclusions
The development, improvement and implementation of fault diagnostics techniques for induction machines requires the use of fast and accurate dynamic models of the machine, which can take into account the asymmetries generated by the machine faults. In this paper, a novel approach has been proposed, using the concept of partial inductance. Instead of using the coil as the basic winding unit, the partial inductance between conductors has been proposed for building the matrices of phase inductances and their derivatives. The partial inductance of a single conductor has been obtained analytically using the magnetic vector potential, and the combination of the partial inductances of all phases has been solved using the FFT, which makes the proposed approach very fast to compute and very easy to implement. The proposed method has been theoretically presented and experimentally validated using the diagnosis of double breakage faults in the squirrel cage of a commercial induction motor for different, non-consecutive positions of the broken bars.
The application of this novel approach is not limited to the specific fault that has been used to perform the experimental validation, namely a broken bars fault. Other types of fault that can be simulated with the procedure presented in this paper include rotor asymmetries of wound-rotor inductions machines, stator inter-turn faults, bearing faults, oscillating loads and supply imbalances. Future work will extend this approach to induction machines with non-uniform air gap, which will allow the simulation of eccentricity faults (static, dynamic and mixed), and the effect of core saturation. In these cases, the analytical computation of phase inductances is a highly complex issue. The proposed method can alleviate this complexity by computing analytically only the MVP generated by a single conductor, which, inserted in the formulae presented in this work, Equations (50) and (51), can provide seamlessly the phases mutual inductances, regardless the complexity of the windings layout. Funding: This work was supported by the Spanish "Ministerio de Economía y Competitividad" in the framework of the "Programa Estatal de Investigación, Desarrollo e Innovación Orientada a los Retos de la Sociedad" (project reference DPI2014-60881-R).

Conflicts of Interest:
The authors declare no conflict of interest.
Machine dimensions: Effective length of the magnetic core = 70.2 mm, radius at the middle of the air gap = 41.1 mm, air gap length = 1.2 mm.