Winding Tensor Approach for the Analytical Computation of the Inductance Matrix in Eccentric Induction Machines

Abstract

Induction machines (IMs) are critical components of many industrial processes, what justifies the use of condition-based maintenance (CBM) systems for detecting their faults at an early stage, in order to avoid costly breakdowns of production lines. The development of CBM systems for IMs relies on the use of fast models that can accurately simulate the machine in faulty conditions. In particular, IM models must be able to reproduce the characteristic harmonics that the IM faults impress in the spatial waves of the air gap magneto-motive force (MMF), due to the complex interactions between spatial and time harmonics. A common type of fault is the eccentricity of the rotor core, which provokes an unbalanced magnetic pull, and can lead to destructive rotor-stator rub. Models developed using the finite element method (FEM) can achieve the required accuracy, but their high computational costs hinder their use in online CBM systems. Analytical models are much faster, but they need an inductance matrix that takes into account the asymmetries generated by the eccentricity fault. Building the inductance matrix for eccentric IMs using traditional techniques, such as the winding function approach (WFA), is a highly complex task, because these functions depend on the combined effect of the winding layout and of the air gap asymmetry. In this paper, a novel method for the fast and simple computation of the inductance matrix for eccentric IMs is presented, which decouples the influence of the air gap asymmetry and of the winding configuration using two independent tensors. It is based on the construction of a primitive inductance tensor, which formulates the eccentricity fault using single conductors as the simplest reference frame; and a winding tensor that converts it into the inductance matrix of a particular machine, taking into account the configuration of the windings. The proposed approach applies routine procedures from tensor algebra for performing such transformation in a simple way. It is theoretically explained and experimentally validated with a commercial induction motor with a mixed eccentricity fault.

1. Introduction

Induction machine (IM) maintenance, integrated in condition-based maintenance (CBM) systems [1,2,3,4,5], is a field of growing industrial interest, due to its widespread use in production lines, electrical vehicles, wind generators, etc. The failure of an IM can cause huge losses, due to unexpected breakdowns of machines and supply systems. To be responsive, CBMs must be able to operate on-line, in a non-invasive way, so that any fault can be detected in an incipient state and corrective measures can be deployed before the fault gets worse [6,7,8,9,10,11]. This requires fast and simple fault diagnostic techniques [12], that can be implemented in embedded field devices, such as digital signal-processors (DSPs) or field-programmable arrays (FPGAs). One of such diagnostic techniques relies on the design of sliding mode observers (SMO) for observing IM states obtained from the healthy and faulty model. In [13], an SMO that uses only input quantity information for on-line broken rotor bar detection is proposed, and in [14], a high-order SMO is designed for detecting inter-turn short circuit faults in IMs. Recent developments [15,16] in this field propose the use of a reduced-order SMO for fault estimation, which is able to simultaneously obtain the exact estimation of state, actuator faults, sensor faults, and extra disturbances. In [16], a novel approach is proposed without employing the equivalent output error injection technology, to overcome the problem of the traditional SMO in application to Markovian jump systems, and in [15], a new proposal is made for avoiding the sliding surface switching problem. Another diagnostic technique is to run an electromechanical model of the machine [17,18,19,20] and compare the simulation outputs (currents and voltages) with the quantities measured at the machine terminals. Divergences between the predicted and the measured values are an indicator of a possible fault, especially if these differences increase over time.

The IM models needed in the aforementioned diagnostic techniques can be built using the finite elements method (FEM) with a very high accuracy [21,22,23], but FEM demands huge computing resources, in terms of time and memory, which hinders its use in low-power embedded units. A faster and leaner alternative is to use analytical models [23,24] that can reproduce the characteristic harmonics induced in the current by a given fault. Another diagnostic area in which IM models are used is in the training of neural networks or expert systems for fault diagnosis [25,26,27,28,29,30,31], which need thousands of tests performed under different working conditions with controlled degrees of a machine fault. In this area, again, the speed of analytical models can give them a decisive advantage over FEM models.

Different IM analytical models for fault diagnosis have been proposed in the technical literature, based on the machine equations expressed in different coordinate systems, such as $d-q$, revolving fields, etc. These models rely on the calculation of the self and mutual inductances between all the machines phases, and their derivatives, as a function of the rotor position. This is a complex, non-linear function, which depends on the windings configurations, and on the rotor position. Besides, in case of a faulty machine, the air gap length or the configuration of the windings may become asymmetrical, making it difficult to use labour-saving procedures that are valid only for symmetrical conditions. In particular, the eccentricity fault [32,33] gives rise to a non-uniform air gap length, which becomes a function of the angular coordinate. Moreover, this function can be different for each rotor position [34]. To overcome this difficulty, the analytical methods for calculating the inductance matrix commonly apply the simplification of considering a sinusoidal distribution of the spatial waves in the air-gap, thus limiting the calculation of the inductances to its fundamental harmonic component. Nevertheless, complex interactions between spatial and time harmonics are present in a faulty machine, but are missing in models restricted to the fundamental component.

Several approaches have been used in the technical literature for obtaining the inductance matrix needed in analytical models. Its components can be determined by direct measurements, as in [35,36], or computed numerically. FEM models have been used for inductance computation in [23], and in [37] a FEM model is combined with a Preisach model for iron loss evaluation. An alternative is to use analytical methods for inductance computation. In [38] a review of the existing methods for the analytical computation of self and mutual inductances in a rotating electrical machine are described, and a new approach based on energy method is presented. A drawback of the analytical methods is that they do not take into account the saturation, the iron path, or the end leakage inductances. In [39] these factors have been simulated via modified air gap length functions.

Instead of a direct, analytical computation of the inductances of phases with a complex winding layout, a successful approach is to start with the inductances of elementary coils, and to combine them via connection matrices to obtain the phase inductances [40]. This approach has been followed also in the winding function approach (WFA) [39]. In [41] the WFA has been combined with a conformal transformation in order to take into account air gap length variations due to the slots. A drawback of these methods is that they need complex winding functions, which depend on the relative position of the coils and the rotor position. Moreover, these functions depend on the combined effect of the winding layout and of the air gap asymmetry, which makes their computation a highly complex task.

In this work, this line of research is followed, with two fundamental novelties: replacing the coils by the conductor as the basic, most simple winding unit, and using routine tensor algebra for the analytical computation of the inductance matrix. This approach allows decoupling the combined effects of the air gap asymmetry and of the winding configuration in the calculation of the inductance matrix, greatly simplifying its analytical computation. The proposed method for calculating the inductance matrix is developed in two steps:

• First, a primitive inductance tensor is calculated in a reference frame that consists of a thin cylindrical sheet of a high number of parallel bars, statically fixed to the air gap. This can be considered, as [42] states, as a canonical coordinate system, in which the components of the primitive inductance tensor are the same for every IM, except for a scaling factor.
• The primitive inductance tensor is transformed into the final one via a winding tensor [40], which contains the current-sheet generated by each phase when fed by a unit current, using routine tensor algebra procedures.

The proposed approach neatly decouples the geometrical configuration of the machine air gap (which can be asymmetrical), represented in the primitive inductance tensor, and the configuration of the phase windings (which can be arbitrarily complex), represented in the winding tensor. Both tensors are defined analytically in an independent way, which simplifies their formulation. Their combined effect, obtained using routine tensor algebra operations, gives the final IM inductance matrix in a simple and fast way, which may be denoted as the winding tensor approach (WTA).

The structure of the paper is the following one. In Section 2, the analytical model of the IM used in this work, in a natural coordinate system, is presented. The analytical computation of the inductance matrix that appears in this model is developed in Section 3 for the case of a healthy and an eccentric IM, using tensor algebra. In Section 4 tensor algebra is applied again to take into account different phase connections, as those imposed by a squirrel cage rotor. An experimental validation of the proposed approach is carried out using a commercial IM motor with a provoked mixed eccentricity fault. This motor is first simulated in Section 5, and the results are compared with those obtained from the experimental tests in Section 6. Finally, Section 7 presents the conclusions of this work.

2. Analytical Model of the IM Using a Natural Coordinate System

Let’s consider a generic IM with $n_s$ stator phases and $n_r$ rotor phases, with a total number of phases $n=n_s+n_r$. Among the multiple coordinate systems that can be used for obtaining an analytical model of the IM ($dq$, symmetrical components, revolving fields, etc.), a natural coordinate system has been chosen in this work: each stator phase ($s_1,s_2\dots s_{n_s}$) has its own axis as a stationary coordinate axis, and each rotor phase ($r_1,r_2\dots r_{n_r}$) has its own axis as a moving coordinate axis, attached to the phase conductors. All the n phase currents in this coordinate system are considered to be independent variables, thus defining an n-dimensional space. From the point of view of the IM simulation, this choice has the advantage of directly giving the phase currents, without needing any further transformation.

In the natural coordinate system, two equations are needed to model the IM operation [40,43]

$$\begin{array}{cccc}•& \mathrm{Equation}\mathrm{of}\mathrm{voltage}:& & \mathit{e}=\mathit{R}\mathit{i}+\mathit{L}\frac{d\mathit{i}}{dt}+\mathit{i}\frac{d\mathit{L}}{d\theta }\dot{\theta }\\ •& \mathrm{Equation}\mathrm{of}\mathrm{torque}:& & T=R_{\theta }\dot{\theta }+J\frac{d\dot{\theta }}{dt}-\frac{1}{2}{\mathit{i}}^t\frac{d\mathit{L}}{d\theta }\mathit{i}\end{array}$$

(1)

where the subscript ‘${}^t$’ stands for the transpose operator. The quantities that appear in (1) are the following ones:

• $\mathit{i}$ is the current tensor. Its components are the instantaneous current in each winding $\mathit{i}={[i_1,i_2,\dots ,i_n]}^t$.
• $\mathit{e}$ is the voltage tensor. Its components are the instantaneous terminal voltages applied to each winding $\mathit{e}={[e_1,e_2,\dots ,e_n]}^t$.
• $\mathit{R}$ is the resistance tensor. Its components are the resistances of all windings. It is a symmetrical dyadic tensor, an square array of $n^2$ constant components.
• $\mathit{L}$ is the inductance tensor. Its components are the self and mutual inductances of all windings along the electrical axes. It is a symmetrical dyadic tensor, an square array of $n^2$ elements. It can be expressed as the sum of two components, one with the inductances corresponding to the main flux linkages ${\mathit{L}}_{\mathit{m}}$, and other with the leakage inductances ${\mathit{L}}_{\mathit{\sigma }}$, as

  $$\mathit{L}={\mathit{L}}_{\mathit{m}}+{\mathit{L}}_{\mathit{\sigma }}$$

  (2)
  End turns, end rings, and slot leakage inductances, included in the ${\mathit{L}}_{\mathit{\sigma }}$ matrix, need to be pre-calculated, as usual in the technical literature, where explicit expressions for these inductances can be found in [44,45,46]. This work deals only with the analytical computation of ${\mathit{L}}_{\mathit{m}}$ in (2). Linear behavior of the iron material will be assumed, as in [47]. This limitation of the analytical model can be overcome using a modified air gap length function to take into the saturation, as in [39].
• The rest of the terms that appear in (1) are the instantaneously applied shaft torque T, the frictional resistance of the shaft $R_{\theta }$, and the moment of inertia J.

The equation of voltage in (1) can be expressed in a more condensed form making use of the tensor of flux linkages of the IM phases, $\mathit{\phi }=\mathit{L}\mathit{i}$. Besides, neglecting the frictional resistance of the shaft ($R_{\theta }$ = 0) in the equation of torque in (1), the set of electro-mechanical equations of the IM, in the natural coordinate system, is given by

$$\left\{\begin{array}{ccc}\mathit{e}& =& \mathit{R}\mathit{i}+\frac{d\mathit{\phi }}{dt}\\ T& =& J\frac{d\dot{\theta }}{dt}-\frac{1}{2}{\mathit{i}}^t\frac{d\mathit{L}}{d\theta }\mathit{i}\end{array}\right.$$

(3)

An implementation of (3) in a Simulink model is shown in Figure 1. The mutual inductances between the stator and rotor phases in ${\mathit{L}}_{\mathit{m}}$ (2) depend on the rotor position, and must be updated at each step of the simulation.

Transformation of the Coordinate System

The quantities $\mathit{i}$, $\mathit{e}$, $\mathit{R}$ and $\mathit{L}$ are tensors, that is, if a different coordinate system is chosen (for example, the hypothetical axis of symmetrical components, or a stationary $dq$ coordinate system), these quantities remain invariant. Only their components (currents, voltages, self and mutual inductances, etc.) in the new coordinate system will be transformed, in the same way that an invariant vector can have different components under different coordinate systems, in spite of not changing neither its modulus nor its orientation. In particular, if the current tensor is expressed in a coordinate system different than the natural one, its new components ${\mathit{i}}^{\prime }$ would be different than the old ones, $\mathit{i}$. Nevertheless, if the matrix $\mathit{C}$ of the coordinate transformation is given, then the relation between the old components and the new ones can be expressed as

$$\mathit{i}=\mathit{C}{\mathit{i}}^{\prime }$$

(4)

and the transformation law of the rest of tensors $\mathit{e}$, $\mathit{R}$ and $\mathit{L}$ is given, applying tensor algebra, by

$$\begin{array}{ccc}{\mathit{e}}^{\prime }& =& {\mathit{C}}^t\mathit{e}\\ {\mathit{L}}^{\prime }& =& {\mathit{C}}^t\mathit{L}\mathit{C}\\ {\mathit{R}}^{\prime }& =& {\mathit{C}}^t\mathit{R}\mathit{C}\end{array}$$

(5)

where ${\mathit{C}}^t$ stands for the transpose of matrix $\mathit{C}$.

3. Computation of the Mutual Inductance Matrix Using Tensor Algebra

Neglecting the iron saturation and losses, mutual inductances depend only on the geometry of the system [48]. Therefore, their computation is done in this work in the spatial domain of the air gap, using the current-sheet generated by the phase currents $\mathit{i}$. The steps of this approach are:

• Definition of the canonical coordinate system for representing the current-sheet distribution along the air gap periphery. In this system, the components of the tensors $\mathit{i}$ and ${\mathit{L}}_{\mathit{m}}$ are independent of the connections of the phase conductors.
• Calculation of the current-sheet from the phase currents $\mathit{i}$. The winding tensor contains the connections between the conductors of the phases, which can be arbitrarily complex.
• Definition of the primitive inductance tensor in the canonical coordinate system, which is independent of the layout of the winding, and the same for every IM, apart from a scaling factor.
• Transformation of the inductance tensor to the natural coordinate system using tensor algebra (5).

3.1. The Components of the Current Tensor

The components of the current tensor are obtained first in the canonical coordinate system, where they are independent of the connections of the phase conductors, and after transformed to the natural coordinate system, using a winding tensor that represents the windings layout.

3.1.1. The Current Tensor in the Canonical Coordinate System

The physical representation of the current tensor $\mathit{i}$ in rotating electrical machines is a current-sheet distributed along the air gap periphery [43]. The most suitable coordinate system to represent it consists in a thin cylindrical sheet of N parallel bars [40], statically placed at the air gap, as shown in Figure 2. The value of N must be high to achieve a high spatial resolution (N = 3600 in this work). The width of each individual bar is assumed to be $2\pi /N$, while its height is considered negligible. An electrical coordinate axis is attached to each conductor, thus defining a N dimensional space in which any current-sheet can be represented with up to $N/2$ spatial harmonics. This N dimensional space is spanned by a basis with N elements, given by

$$\left\{\begin{array}{ccc}{\mathit{u}}_{\mathit{c}\mathbf{1}}& =& {[1,0,0,\dots ,0]}^t\\ {\mathit{u}}_{\mathit{c}\mathbf{2}}& =& {[0,1,0,\dots ,0]}^t\\ & \dots & \\ {\mathit{u}}_{\mathit{cN}}& =& {[0,0,0,\dots ,1]}^t\end{array}\right.$$

(6)

where the kth basis element ${\mathit{u}}_{\mathit{ck}}$ has all components equal to 0 except the kth that is 1. The vectors ${\mathit{u}}_{\mathit{ck}}$ are unitary and orthogonal, and they form an ordered basis, which is called the standard or canonical basis. In this basis, the current-sheet can be represented as a linear combination of the basis elements as

$${\mathit{i}}_{\mathit{c}}=\sum _{k=1}^N{i}_{ck}·{\mathit{u}}_{\mathit{ck}}$$

(7)

where $i_{ck}$ the current in conductor k. That is, the components of the current tensor in this coordinate system, ${\mathit{i}}_{\mathit{c}}$, are simply the N currents through the N independent conductors

$${\mathit{i}}_{\mathit{c}}={[i_{c_1},i_{c_2},\dots ,i_{c_N}]}^t$$

(8)

3.1.2. Transformation of the Current Tensor to Natural Coordinates

The current-sheet ${\mathit{i}}_{\mathit{c}}$ (8) can be also expressed in the natural coordinate system, using the n phase currents as independent variables. This n dimensional space is spanned by a basis formed by n basis vectors, $n_s$ stator and $n_r$ rotor vectors. In this basis, the current-sheet ${\mathit{i}}_{\mathit{c}}$ can be represented as a linear combination of the new basis vectors as

$${\mathit{i}}_{\mathit{c}}=\sum _{k=1}^n{i}_k·{\mathit{z}}_{\mathit{k}}$$

(9)

where $i_k$ represents the current in phase k. Each basis vector in (9) has N components, which for the kth basis vector ${\mathit{z}}_{\mathit{k}}$ are the ampere-turns generated by phase k at each angular interval of Figure 2, when fed with a unit dc current. This value coincides with the number of conductors of phase k in each interval, with a ± sign corresponding to the direction of the current.

$$\left\{\begin{array}{ccc}{\mathit{z}}_{\mathit{1}}& =& {[z_{11},z_{21},z_{31},\dots ,z_{N1}]}^t\hfill \\ {\mathit{z}}_{\mathit{2}}& =& {[z_{12},z_{22},z_{32},\dots ,z_{N2}]}^t\hfill \\ & ⋮& \\ {\mathit{z}}_{\mathit{n}}& =& {[z_{1n},z_{2n},z_{3n},\dots ,z_{Nn}]}^t\hfill \end{array}\right.$$

(10)

The number of basis vectors in this new coordinate system (10) is much lower than in the primitive coordinate system (6). Nevertheless, they are neither unitary nor orthogonal.

The new vector basis (10) can be expressed in the canonical basis (6) as

$${\mathit{z}}_{\mathit{k}}=\sum _{l=1}^N{z}_{kl}·{\mathit{u}}_{\mathit{l}}$$

(11)

Using (10) and (11), (9) becomes

$${\mathit{i}}_{\mathit{c}}=\sum _{k=1}^n\sum _{l=1}^N{i}_k·z_{kl}·{\mathit{u}}_{\mathit{l}}$$

(12)

This coordinate transformation can be formulated using a $(N\times n)$ transformation matrix ${\mathit{C}}_{\mathit{c}}$ as (4)

$${\mathit{i}}_{\mathit{c}}={\mathit{C}}_{\mathit{c}}\mathit{i}$$

(13)

where the columns of ${\mathit{C}}_{\mathit{c}}$ are the new basis vectors (10),

$${\mathit{C}}_{\mathit{c}}=\left[\begin{array}{cccc}{z}_{11}& z_{12}& \cdots & z_{1n}\\ z_{21}& z_{22}& \cdots & z_{2n}\\ ⋮& ⋮& \ddots & ⋮\\ z_{N1}& z_{N2}& \cdots & z_{Nn}\end{array}\right]$$

(14)

The transformation matrix ${\mathit{C}}_{\mathit{c}}$ represents the current constraints imposed by the connections between the conductors of each winding. Therefore, as Kron states in [40], this particular transformation matrix can be considered as one aspect of the transformation tensor ${\mathit{C}}_{\mathit{c}}$, that will be referred to as the winding tensor. Its $(i,j)$ element contains the number of phase conductors of phase j in an angular interval $\pi /2N$, centered at $i·\frac{2\pi }{N}$. This winding tensor must be obtained for the N possible angular positions of the rotor (${\theta }_k=k·\frac{2\pi }{N}$, with $k=0,\dots ,N-1$). Nevertheless, the columns of ${\mathit{C}}_{\mathit{c}}$ corresponding to the rotor phases for a given rotor position ${\theta }_k$ are the same as the columns defined with the rotor at the origin (${\theta }_0=0$), but rotated k positions.

3.1.3. The Winding Tensor for Phases with the Same Configuration

In (14) no restrictions are imposed on the connections of the conductors of each phase, which can be arbitrarily complex, as in the case of asymmetrical windings (turn-to-turn short circuits, broken bars, etc.). Nevertheless, in case of a healthy machine, the configuration of all the phases of a particular winding (stator or rotor) is the same. Therefore, the vector column of ${\mathit{C}}_{\mathit{c}}$ corresponding to the kth stator phase is equal to the vector column of the first stator phase, but rotated $k·N/n_s$ elements to the bottom. The same applies to the rotor phases, but in this case the rotation is $k·N/n_r$ positions.

3.2. The Components of the Mutual Inductance Tensor

The components of the mutual inductance tensor will be obtained first in the canonical coordinate system, where they are independent of the connections of the phase conductors, and then transformed to the natural coordinate system, using the winding tensor.

3.2.1. The Primitive Inductance Tensor

In the canonical coordinate system (6), the mutual inductance tensor, ${\mathit{L}}_{\mathit{mc}}$, is a $N\times N$ square matrix

$${\mathit{L}}_{\mathit{mc}}=\left[\begin{array}{cccc}{L}_{m{c}_{11}}& L_{m{c}_{12}}& \cdots & L_{m{c}_{1N}}\\ L_{m{c}_{21}}& L_{m{c}_{22}}& \cdots & L_{m{c}_{2N}}\\ ⋮& ⋮& \ddots & ⋮\\ L_{m{c}_{N1}}& L_{m{c}_{N2}}& \cdots & L_{m{c}_{NN}}\end{array}\right]$$

(15)

whose component $(i,j)$, $L_{m{c}_{ij}}$, is the mutual partial inductance [24] between the conductors placed at positions $i·\frac{2\pi }{N}$ and $j·\frac{2\pi }{N}$. The tensor ${\mathit{L}}_{\mathit{mc}}$ (15) will be denoted as the primitive inductance tensor.

3.2.2. The Primitive Inductance Tensor of a Non-Eccentric IM

In case of IMs with uniform air gap, as represented in Figure 2, $L_{m{c}_{ij}}$ depends only on the angular separation between conductors i and j, as given in [49]

$${\mathit{L}}_{\mathit{mc}}(i,j)=L_{m{c}_{ij}}=\frac{{\mu }_0·l·r·\pi }{g}·{\left(\frac{1}{2}-\frac{|i-j|}{N}\right)}^2$$

(16)

where ${\mu }_0=4\pi \times {10}^{-7}$ T·m·A⁻¹, l is the effective length of the stator bore, r is the radius at the center of the air gap, and g is the air gap length.

From (16), the components of ${\mathit{L}}_{\mathit{mc}}$ are the same for every IM, except for the scaling factor $\frac{{\mu }_0·l·r·\pi }{g}$, which depends only on the geometrical dimensions of the machine, l, r and g. Besides, ${\mathit{L}}_{\mathit{mc}}$ is a circulant, symmetrical matrix, where every column vector is obtained by rotating one element to the bottom of the preceding column vector.

3.2.3. The Primitive Inductance Tensor of an Eccentric IM

In the case of rotor eccentricity the air gap length is not uniform, because the rotor center $O_r$ does not coincide with the stator center $O_s$, as shown in Figure 3.

From Figure 3, the position of the rotor center can be represented using its radial coordinate, ${\delta }_r·g_0$, and its angular coordinate ${\mathsf{\Theta }}_r$, as

$$\overrightarrow{O_s{O}_r}=g_0·{\delta }_r·e^{j{\mathsf{\Theta }}_r}0\le {\delta }_r<1,0\le {\mathsf{\Theta }}_r<2\pi$$

(17)

where $g_0$ is the air gap length of the IM without any eccentricity, and ${\delta }_r$ is the degree of eccentricity ($0\le {\delta }_r<1$). Additionally, the coordinates ($g_0·{\delta }_r$, ${\mathsf{\Theta }}_r$) of the rotor center depend on the angular position of the rotor ${\theta }_r$, and the degree of static ${\delta }_{se}$ and dynamic ${\delta }_{de}$ eccentricity of the machine (see Figure 4), as

$${\mathsf{\Theta }}_r\left({\theta }_r\right)={tan}^{-1}\left(\frac{{\delta }_{de}sin\left({\theta }_r\right)}{{\delta }_{se}+{\delta }_{de}cos\left({\theta }_r\right)}\right)$$

(18)

$${\delta }_r\left({\theta }_r\right)=\sqrt{{{\delta }_{se}}^2+{{\delta }_{de}}^2+2{\delta }_{se}{\delta }_{de}cos\left({\theta }_r\right)}$$

(19)

where ${\theta }_r$ represents the angle of rotation of the machine rotor.

For computing the inductance matrix, the inverse of the air gap length function is needed to obtain the permeance function of the machine. It can be fully defined in terms of the coordinates of the rotor center (17) as [23]

$$g{(\phi ,{\mathsf{\Theta }}_r,{\delta }_r)}^{-1}=g_0^{-1}·\left(A_0+\sum _{m=1}^{n_t}{A}_m·cos\left(m(\phi -{\mathsf{\Theta }}_r)\right)\right)$$

(20)

where

$$A_0=\frac{1}{\sqrt{1-{{\delta }_r}^2}}{A}_m=2{\left(\frac{1-\sqrt{1-{{\delta }_r}^2}}{\sqrt{1-{{\delta }_r}^2}}\right)}^{m}m=1\dots n_t$$

(21)

It is worth mentioning that only the first term of the series in (20) has been used in [50,51,52,53], and two terms in [54]. In this paper, (20) can take into account a generic number $n_t$ of terms, where the value of $n_t$ can be chosen to achieve the desired precision.

Each component (i,j) of the induction matrix $L_{m{c}_{ij}}$ in an eccentric IM depends not only on the angular separation between conductors i and j, but also on their absolute position and on the position of the rotor center, whose coordinates ($g_0·{\delta }_r$, ${\mathsf{\Theta }}_r$) are, in turn, functions of the rotor angular position (18), (19). The analytical expression of $L_{m{c}_{ij}}$, for a given rotor position ${\theta }_k=k·\frac{2\pi }{N}$, can be expressed as [24]

$${\mathit{L}}_{\mathit{mc}}(i,j){|}_k=\frac{{\mu }_{0}lr}{g_0}·\mathsf{\Lambda }(i\frac{2\pi }{N},j\frac{2\pi }{N},{\mathsf{\Theta }}_r\left(k\frac{2\pi }{N}\right),{\delta }_r\left(k\frac{2\pi }{N}\right))$$

(22)

where

$$\begin{array}{c}\mathsf{\Lambda }(\alpha ,\phi ,{\mathsf{\Theta }}_r,{\delta }_r)=\frac{A_0}{4\pi }{(\phi -\alpha )}^2+{\displaystyle \sum _{m=1}^{n_t}}\frac{A_m}{2\pi }\left(\frac{(\phi -\alpha )sin\left(m(\phi -{\mathsf{\Theta }}_r)\right)}{m}+\frac{cos\left(m(\phi -{\mathsf{\Theta }}_r)\right)}{m^2}\right)-\\ \left(\frac{1}{2}-K(\alpha ,{\mathsf{\Theta }}_r,{\delta }_r)\right)·\left(A_0(\phi -\alpha )+{\displaystyle \sum _{m=1}^{n_t}}{A}_m\frac{sin\left(m(\phi -{\mathsf{\Theta }}_r)\right)}{m}\right)\end{array}$$

(23)

and

$$K(\alpha ,{\mathsf{\Theta }}_r,{\delta }_r)=\sum _{m=1}^{n_t}\frac{A_m}{2\pi A_0}\frac{sin\left(m({\mathsf{\Theta }}_r-\alpha )\right)}{m}$$

(24)

From (22), the components of ${\mathit{L}}_{\mathit{mc}}$ are the same for every IM with a given degree of static and dynamic eccentricity, except for the scaling factor $\frac{{\mu }_0·l·r·\pi }{g_0}$, which depends only on the geometrical dimensions of the machine.

It is worth remarking that the primitive inductance tensor ${\mathit{L}}_{\mathit{mc}}$ includes the effect of the air gap asymmetry generated by the mixed eccentricity fault, but is independent of the winding configuration, because it has been obtained using the conductor as the basic unit. Therefore, it is valid for any IM, except for a scaling factor. This leads to a great simplification compared with other existing methods, such as the WFA, which rely on winding functions whose definition depends both on the air gap asymmetry and on the configuration of the winding coils.

3.2.4. Transformation of the Inductance Tensor to Natural Coordinates

In the natural coordinate system, for each rotor position, the inductance tensor ${\mathit{L}}_{\mathit{m}}$ in (2) can be obtained from the primitive inductance tensor ${\mathit{L}}_{\mathit{mc}}$, either considering healthy (16) or eccentric machine (22), and from the winding tensor ${\mathit{C}}_{\mathit{c}}$, as (5)

$${\mathit{L}}_{\mathit{m}}={{\mathit{C}}_{\mathit{c}}}^t{\mathit{L}}_{\mathit{mc}}{\mathit{C}}_{\mathit{c}}$$

(25)

Again, at this point it is worth remarking that the transformation from the primitive inductance tensor ${\mathit{L}}_{\mathit{mc}}$, which includes only the effect of the air gap asymmetry generated by the mixed eccentricity fault, into the final inductance matrix ${\mathit{L}}_{\mathit{m}}$, taking into account the configuration of the windings, is a routine tensor algebra operation (25) that simply consists in multiplying the primitive inductance tensor by the winding tensor ${\mathit{C}}_{\mathit{c}}$ (14) (and its transpose). The winding tensor has been defined without any relation to the air gap asymmetry, simply indicating the number and direction of the winding conductors at each interval of the rotor and stator periphery.

This advantage can be further exploited to introduce the effect of winding related faults in the model through the winding tensor ${\mathit{C}}_{\mathit{c}}$, such as inter-turn short circuits or phase asymmetries. This would allow the analysis of combined eccentricity and winding faults with a small increase in analytical and computational complexity. Due to space constraints, this approach has not been considered in this work, and will be presented in a future one.

4. Additional Current Constraints Imposed by Phase Connections

In (3) the phase currents $\mathit{i}$ are independent variables. Nevertheless, the connections between the phases can make some currents dependent on others. These constraints can be expressed as a connection tensor ${\mathit{C}}_{\mathit{i}}$, which relates the original, independent phase currents $\mathit{i}$ (before considering the interconnections between the phases), with the new ones ${\mathit{i}}^{\prime }$, including the constraints introduced by the phase connections as

$$\mathit{i}={\mathit{C}}_{\mathit{i}}{\mathit{i}}^{\prime }$$

(26)

As the transformation given by ${\mathit{C}}_{\mathit{i}}$ is holonomic, because its components do not depend on the rotor position, (1) remains valid when expressed in terms of the new, reduced quantities ${\mathit{i}}^{\prime }$, ${\mathit{e}}^{\prime }$, ${\mathit{L}}^{\prime }$ and ${\mathit{R}}^{\prime }$, giving

$$\begin{array}{c}\left\{\begin{array}{ccc}{\mathit{e}}^{\prime }& =& {\mathit{R}}^{\prime }{\mathit{i}}^{\prime }+{\mathit{L}}^{\prime }\frac{d{\mathit{i}}^{\prime }}{dt}+{\mathit{i}}^{\prime }\frac{d{\mathit{L}}^{\prime }}{d\theta }\dot{\theta }\\ T& =& R_{\theta }\dot{\theta }+J\frac{d\dot{\theta }}{dt}-\frac{1}{2}{\mathit{i}}^{\prime t}\frac{d{\mathit{L}}^{\prime }}{d\theta }{\mathit{i}}^{\prime }\end{array}\right.\end{array}$$

(27)

where ${\mathit{e}}^{\prime }={{\mathit{C}}_{\mathit{i}}}^t\mathit{e}$, ${\mathit{R}}^{\prime }={{\mathit{C}}_{\mathit{i}}}^t\mathit{R}{\mathit{C}}_{\mathit{i}}$, and ${\mathit{L}}^{\prime }={{\mathit{C}}_{\mathit{i}}}^t\mathit{L}{\mathit{C}}_{\mathit{i}}$.

A particular example of current constraints in squirrel cage IMs are those imposed by the physical configuration of the rotor cage, which are analyzed in the following Section.

Current Constraints in a Squirrel Cage IM

In the analytical model of a squirrel cage IM, the $n_s$ stator phases are considered as independent electrical circuits in (3), without any current constraints. Therefore, the resistance and leakage tensors of the stator winding are square matrices of size $n_s\times n_s$, whose diagonal terms are equal to the stator phase resistances, $R_s$, and to the phase leakage inductances, $L_{\sigma s}$, respectively. All the terms outside the diagonals are zero. If all the stator phases have the same configuration then

$${\mathit{R}}_{\mathit{s}}=\left[\begin{array}{ccc}{R}_s& & 0\\ & \ddots & \\ 0& & R_s\end{array}\right]$$

(28)

$${\mathit{L}}_{\mathit{\sigma }\mathit{s}}=\left[\begin{array}{ccc}{L}_{\sigma s}& & 0\\ & \ddots & \\ 0& & L_{\sigma s}\end{array}\right]$$

(29)

where $R_s$ and $L_{\sigma s}$ are the resistance and the leakage inductance of a stator phase, respectively.

The electric circuit of the squirrel cage rotor, with $n_b$ bars, can be built using $n_b$ rotor loops (each loop formed by two consecutive rotor bars), plus two additional loops corresponding to the end rings [55,56], as represented in Figure 5. The currents in the loops formed by two consecutive bars ($i_{b1}$ to $i_{b{n}_b}$) are coupled to each other and to the stator currents through their mutual inductances. On the contrary, the end ring loop currents ($i_{e1}$ and $i_{e2}$ in Figure 5) do not couple with the stator currents, and couple with the other rotor loop currents only through the end ring leakage inductances and the end ring resistances ($L_{\sigma e}$ and $R_e$ in Figure 5, respectively). Therefore, the resistance and leakage tensors of the squirrel cage rotor (with $n_b$ bars) are square matrices of size $(n_b+2)\times (n_b+2)$, due to the presence of the two extra loops formed by the end rings.

The resistance matrix of the rotor cage of Figure 5 is given by

$${\mathit{R}}_{\mathit{b}}=\left[\begin{array}{ccccccccc}{R}_{be}& -R_b& & & & & & -R_e& -R_e\\ -R_b& R_{be}& -R_b& & & & & -R_e& -R_e\\ & -R_b& R_{be}& -R_b& & 0& & -R_e& -R_e\\ & & \ddots & \ddots & \ddots & & & ⋮& ⋮\\ & 0& & -R_b& R_{be}& -R_b& & -R_e& -R_e\\ & & & & -R_b& R_{be}& -R_b& -R_e& -R_e\\ & & & & & -R_b& R_{be}& -R_e& -R_e\\ -R_e& -R_e& \dots & -R_e& -R_e& -R_e& -R_e& n_b{R}_e& 0\\ -R_e& -R_e& \dots & -R_e& -R_e& -R_e& -R_e& 0& n_b{R}_e\end{array}\right]$$

(30)

with $R_{be}=2(R_b+R_e)$.

Nevertheless, it is advisable to reduce the set of rotor currents to those circulating in the loops that contain rotor bars, which are the only ones coupled through the mutual loop inductances ${\mathit{L}}_{\mathit{m}}$ (25). In this way, the tensor of rotor currents is reduced to

$${\mathit{i}}_{\mathit{r}}={[i_{b1},i_{b2},\dots ,i_{b{n}_b-1}]}^t$$

(31)

This reduction can be achieved because there are two current constraints in the electrical circuit Figure 5. In effect, if the end ring loops mutual inductances with the rest of the rotor loops and with the stator phases are neglected, being mainly end winding flux linkages, the currents in the end ring loops can be expressed as:

$$n_b\left(R_e{i}_{e1}+L_{\sigma e}\frac{d{i}_{e1}}{dt}\right)-\sum _{i=1}^{n_b-1}(R_e{i}_{bi}+L_{\sigma e}\frac{d{i}_{bi}}{dt})=0\Rightarrow i_{e1}=\frac{1}{n_b}\sum _{i=1}^{n_b-1}{i}_{bi}$$

(32)

$$n_b\left(R_e{i}_{e2}+L_{\sigma e}\frac{d{i}_{e2}}{dt}\right)-\sum _{i=1}^{n_b-1}(R_e{i}_{bi}+L_{\sigma e}\frac{d{i}_{bi}}{dt})=0\Rightarrow i_{e2}=\frac{1}{n_b}\sum _{i=1}^{n_b-1}{i}_{bi}$$

(33)

These constraints can be formulated as in (26). The old set ${\mathit{i}}_{\mathit{b}}$ of $n_b+1$ rotor currents ($n_b-1$ current loops that include the bars plus two end ring currents) can be obtained from the new set ${\mathit{i}}_{\mathit{r}}$ of $n_b-1$ currents (just the current loops containing the bars), using (32) and (33), as

$$\begin{array}{ccccc}\underbrace{\left[\begin{array}{c}{i}_{b1}\\ i_{b2}\\ ⋮\\ i_{b{n}_b-1}\\ i_{e_1}\\ i_{e_2}\end{array}\right]}& =& \underbrace{\left[\begin{array}{cccc}1& & & \\ & 1& 0& \\ & 0& \ddots & \\ & & & 1\\ \frac{1}{n_b}& \frac{1}{n_b}& \dots & \frac{1}{n_b}\\ \frac{1}{n_b}& \frac{1}{n_b}& \dots & \frac{1}{n_b}\end{array}\right]}& ·& \underbrace{\left[\begin{array}{c}{i}_{b1}\\ i_{b2}\\ ⋮\\ i_{b{n}_b-1}\end{array}\right]}\\ {\mathit{i}}_{\mathit{b}}& =& {\mathit{C}}_{\mathit{i}}& ·& {\mathit{i}}_{\mathit{r}}\\ (n_b+1)\times 1& & (n_b+1)\times (n_b-1)& & (n_b-1)\times 1\end{array}$$

(34)

Using the circuit transformation tensor ${\mathit{C}}_{\mathit{i}}$ in (34), which transforms branch resistances into loop ones, as in [57,58], the resistance matrix of the rotor ${\mathit{R}}_{\mathit{r}}$ that includes the current constraints (32) and (33) is

$${\mathit{R}}_{\mathit{r}}={{\mathit{C}}_{\mathit{i}}}^t·{\mathit{R}}_{\mathit{b}}·{\mathit{C}}_{\mathit{i}}$$

(35)

giving the final result

$${\mathit{R}}_{\mathit{r}}=\left[\begin{array}{ccccccc}{R}_{be}& -R_b& & & & & \\ -R_b& R_{be}& -R_b& & & & \\ & -R_b& R_{be}& -R_b& & 0& \\ & & \ddots & \ddots & \ddots & & \\ & 0& & -R_b& R_{be}& -R_b\\ & & & & -R_b& R_{be}& -R_b\\ & & & & & -R_b& R_{be}\end{array}\right]-\frac{2R_e}{n_b}·\left[\begin{array}{cccc}1& 1& \dots & 1\\ 1& 1& \dots & 1\\ ⋮& ⋮& \ddots & ⋮\\ 1& 1& \dots & 1\end{array}\right]$$

(36)

The term $\frac{2R_e}{n_b}$ is small for a high value of $n_b$, which may justify considering only one end ring current loop, as in [59], or even to neglect both end ring current loops.

In a similar way, the matrix of rotor leakage inductances is given by

$$\begin{array}{cc}\hfill {\mathit{L}}_{\mathit{\sigma }\mathit{r}}& =\left[\begin{array}{ccccccc}{L}_{\sigma be}& -L_{\sigma b}& & & & & \\ -L_{\sigma b}& L_{\sigma be}& -L_{\sigma b}& & & & \\ & -L_{\sigma b}& L_{\sigma be}& -L_{\sigma b}& & 0& \\ & & \ddots & \ddots & \ddots & & \\ & 0& & -L_{\sigma b}& L_{\sigma be}& -L_{\sigma b}& \\ & & & & -L_{\sigma b}& L_{\sigma be}& -L_{\sigma b}\\ & & & & & -L_{\sigma b}& L_{\sigma be}\end{array}\right]-\hfill \\ & -\frac{2L_{\sigma e}}{n_b}·\left[\begin{array}{cccc}1& 1& \dots & 1\\ 1& 1& \dots & 1\\ ⋮& ⋮& \ddots & ⋮\\ 1& 1& \dots & 1\end{array}\right]\hfill \end{array}$$

(37)

where $L_{\sigma be}=2(L_{\sigma b}+L_{\sigma e})$.

Using (28), (29), (36) and (37), the final resistance and leakage inductance matrices that are used in the electromechanical equations of the healthy squirrel cage IM are

$$\mathit{R}=\left[\begin{array}{cc}{\mathit{R}}_{\mathit{s}}& 0\\ 0& {\mathit{R}}_{\mathit{r}}\end{array}\right]$$

(38)

$${\mathit{L}}_{\mathit{\sigma }}=\left[\begin{array}{cc}{\mathit{L}}_{\mathit{\sigma }\mathit{s}}& 0\\ 0& {\mathit{L}}_{\mathit{\sigma }\mathit{r}}\end{array}\right]$$

(39)

5. Analytical Model of the Tested IM

In this Section, the analytical model of the commercial IM whose characteristics are given in Appendix A, which is used for the experimental validation of the proposed approach, is calculated considering two different motor conditions:

• Healthy conditions.
• Faulty conditions, with a mixed eccentricity fault, with 30% of static eccentricity (${\delta }_{se}$ = 0.3) and 30% of dynamic eccentricity (${\delta }_{de}$ = 0.3).

For building the analytical model of this motor, first the number N of intervals in which the air gap periphery is divided must be selected. In this work, it has been chosen $N=3600$, giving an angular resolution of ${0.1}^{\circ }$.

5.1. Analytical Model of the Tested IM in Healthy Conditions

5.1.1. Primitive Inductance Tensor of the Healthy IM

The primitive inductance tensor (15) has been computed for the healthy motor using (16), and its first column is represented in Figure 6. Each one of the rest of the columns of (15) is equal to the previous one, with its elements rotated by one position. This tensor is independent of the angular position of the rotor.

5.1.2. Winding Tensor of the Healthy IM

The winding tensor, which contains the current-sheet generated by each phase, when fed with a unit current, has been obtained using the data provided in Appendix A. Figure 7 shows the three first columns of the winding tensor ${\mathit{C}}_{\mathit{c}}$ (14), corresponding to the stator phases, and Figure 8 shows the next three columns of ${\mathit{C}}_{\mathit{c}}$ (14), corresponding to the three first rotor loops. The inclination of rotors slots has been taken into account as in [49].

5.1.3. Mutual Inductance Matrix of the Healthy IM

The mutual inductance matrix ${\mathit{L}}_{\mathit{m}}$ is found applying (15), using the primitive inductance tensor (Figure 9) and the winding tensor (Figure 7 and Figure 8). Figure 9 represents, as a function of the rotor position, the self inductance of the first stator phase (Figure 9, top), of the first rotor loop (Figure 9, middle), and the mutual inductance between the first stator phase and the first rotor loop (Figure 9, bottom). Figure 10 shows their respective angular derivatives.

5.1.4. Resistance and Leakage Inductance Matrices of the Healthy Machine

The resistance (38) and leakage inductance (39) matrices of the healthy machine are assembled using (28), (29), (36) and (37), using the values of the resistance and the leakage inductance of a stator phase, a rotor bar and end ring segment given in Appendix A.

5.2. Analytical Model of the Tested IM with an Eccentricity Fault

The model of the motor given in Appendix A is obtained in this section assuming a mixed eccentricity fault, with a level of 30% of static eccentricity and 30% of dynamic eccentricity (${\delta }_{se}$ = 0.3, ${\delta }_{de}$ = 0.3). It is modeled using the expression for the primitive inductance tensor of the eccentric IM (22).

5.2.1. Primitive Inductance Tensor of the Eccentric IM

The primitive inductance tensor (15) of the eccentric IM has been computed using (22), using ${\delta }_{se}$ = 0.3 and ${\delta }_{de}$ = 0.3 in (18) and (19). Contrary to the case of the healthy machine, the primitive inductance tensor depends now on the rotor position. The columns of this tensor cannot be obtained by a rotation of the previous ones, due to the eccentricity of the rotor. Figure 11 shows the first column of the primitive inductance tensor, for a rotor angular position equal to zero. It is worth mentioning that the primitive inductance tensor just captures the eccentricity fault, being independent of the winding configuration.

5.2.2. Winding Tensor of the Eccentric IM

The winding tensor of the healthy IM, shown in Figure 7 and Figure 8, is not affected by the eccentricity fault. Therefore, it is the same as the winding tensor of the healthy machine.

5.2.3. Mutual Inductance Matrix of the Eccentric IM

The inductance matrix of the faulty IM is found by applying (15), using the primitive inductance tensor represented in Figure 12 and the winding tensor represented in Figure 7 and Figure 8. Figure 12 represents, as a function of the rotor position, the self inductance of the first stator phase (Figure 12, top), of the first rotor loop (Figure 12, middle), and the mutual inductance between the first stator phase and the first rotor loop (Figure 12, bottom). Figure 13 shows their respective angular derivatives.

5.2.4. Resistance and Leakage Inductance Matrices of the Eccentric Machine

The resistance (38) and leakage inductance (39) matrices of the healthy machine are not affected by the eccentricity fault.

6. Experimental Validation

The experimental procedure carried out for the validation of the proposed method consists in provoking an artificial eccentricity fault to the tested motor (Appendix A), and comparing the characteristic fault harmonics that this fault induces in the motor current with those obtained from the simulated motor. To this end, the original bearings of the motor (see Figure 14a) have been replaced with new bearings (Figure 14d) having a smaller outer diameter and a greater inner diameter. These new bearings have been displaced from the center of the stator bore using two precision eccentric steel rings (Figure 14b,c), placed in the bearings housing (Figure 14b) and on the shaft (Figure 14c). The cylindrical surfaces of both rings are eccentric, 0.09 mm in the case of the outer ring b, and 0.09 mm in the case of the inner ring c. This assembly (Figure 14e) results in a rotor with a 30% of static eccentricity and a 30% of dynamic eccentricity.

A mixed eccentricity fault [51] induces two characteristic series of harmonic components in the motor current spectrum: one as side bands around the principal slot harmonics, and other one around the fundamental component. The frequencies of this low frequency series depend on the rotor speed, and can be obtained as

$$f_{ME}\left(s\right)=f_1\pm (k(1-s)f_1/p),k=1,2,3\dots$$

(40)

where $f_1$ is the power supply frequency, s is the slip and p is the number of pole pairs of the machine.

Using the dominant component of the series (40) (k = 1), the mixed eccentricity fault can be detected through the presence in the stator current spectrum of harmonic components at frequencies:

$$f_{ME}\left(s\right)=f_1\pm (1-s)f_1/p=f_1\pm f_r,$$

(41)

where $f_r$ is the rotational frequency of the motor. For the tested motor ($p=2$), this gives

$$f_{ME}\left(s\right)=f_1\pm (1-s)f_1/2$$

(42)

To verify the validity of the method proposed in this paper, in particular its ability to reproduce the fault harmonics at frequencies given by (42), the eccentric motor has been tested during a start-up transient, with a final permanent regime speed of 1445 rpm ($s=(1500-1445)/1500=0.0367$). To this end, one of the phase currents has been sampled, using the current clamp whose data is given in Appendix B, during an acquisition time of 8 seconds, with a sampling rate of 2 kHz. The spectrogram of this current, obtained with the computer platform given in Appendix C, is shown in Figure 15. As given by (42), two fault-related harmonics appear in permanent regime at frequencies $f_{ME}\left(0.0367\right)=50\pm (1-0.0367)·50/2=$ [25.92 Hz, 74.08 Hz]. Both harmonic components have the same frequency as the fundamental component at the beginning of the start-up transient, and their frequency varies proportionally to the speed up to their final frequency in the permanent regime. For comparison purposes, Figure 16 presents the spectrogram obtained during the start-up transient of the same motor in healthy conditions, prior to provoking the eccentricity fault. In this spectrogram there is no presence of the harmonic components produced by the mixed eccentricity fault.

The motor given in Appendix A has been simulated under the same conditions as the experimental test, using the Simulink model given in Figure 1. The spectrogram of the simulated phase current, shown in Figure 17, displays correctly the characteristic signature of the eccentricity fault harmonics in the time-frequency plane, which assesses the validity of the method presented in this work. For comparison purposes, the motor has been simulated in healthy conditions, and in the spectrogram of its current, shown in Figure 18, there is no presence of the harmonic components produced by the mixed eccentricity fault.

7. Conclusions

Tensor algebra provides a powerful tool for develo** the analytical model of IMs, because it allows us to gradually introduce the effects of any fault first at the conductor level, using the primitive inductance tensor, and after at the winding level, using the winding tensor. Routine procedures of tensor algebra facilitate the conversion of the primitive inductance tensor into the inductance matrix of the machine, for any winding configuration. In this work, this new method has been presented, applied to the development of the analytical model of an eccentric IM, and validated using experimental tests. The application of the winding tensor approach for building the analytical model of IMs with other types of faults is currently a work in progress.