On the Mathematical Modeling of Line-Start Permanent Magnet Synchronous Motors under Static Eccentricity

Abstract

Line start permanent magnet synchronous motors experience different types of failures, including static eccentricity. The first step in detecting such failures is the mathematical modeling of the motor under healthy and failed conditions. In this paper, an attempt to develop an accurate mathematical model for this motor under static eccentricity is presented. The model is based on the modified winding function method and coupled magnetic circuits approach. The model parameters are calculated directly from the motor winding layout and its geometry. Static eccentricity effects are considered in the motor inductances calculation. The performance of the line start permanent magnet synchronous motor using the developed mathematical model is investigated using MATLAB/SIMULINK^® software (2013b, MathWorks, Natick, MA, USA) under healthy and static eccentricity condition for different loading values. A finite element method analysis is conducted to verify the mathematical model results, using the commercial JMAG^® software (16.0.02n, JSOL Corporation, Tokyo, Japan). The results show a fine agreement between JMAG^® and the developed mathematical model simulation results.

1. Introduction

Electrical motors play a major role in industrial sectors around the world. They are the backbone of industry and manufacturing. Statistics show that electrical motors consume about two-thirds of the total industrial power consumption in each society [1]. Induction motors constitute by far the largest portion of electrical motors in the market. However, more efficient motors are beginning to appear as other alternatives. A recently manufactured motor, which is used in industrial practices, is the Line-Start Permanent Magnet Synchronous (LSPMS) motor. This motor has a self-starting capability as induction motors due to rotor cage and attains the synchronous speed due to rotor magnets. Under normal operating conditions, a great deal of research work has investigated the LSPMS motor in terms of modeling, efficiency, and energy consumption [2]. The researcher also presents the challenges that are faced by this motor, such as the breaking torque and construction cost [3]. It was shown that this motor enjoys a higher efficiency than induction motors, which in turn, reflected on the energy and money saving benefits [4,5]. As a result, this will open the door for greater penetration in the market. H. Behbahanifard et al. [6] introduced a comprehensive survey of the work done on the modeling trends of healthy LSPMS motor, which includes the initial start and synchronization process, cogging torque, and the armature reaction during the starting process, which results in permanent magnet demagnetization. Kahourzade et al. [7] focused on modeling with more emphasis on the dynamic and steady state behavior of the motor where the LSPMS motor (3-phase, 4-pole, 2 kW) was compared with a similar rating induction motor. It was found that the LSPMS motor had higher efficiency and can withstand overload conditions more than the induction motors. Besides, the starting current for the LSPMS motor is lower than the induction motor at the same load condition. On the same manner, a high performance LSPMS motor model was developed from an existing induction motor model in [5]. Dan Stoia et al. [8] proposed an analytical method to analyze a LSPMS motor based on a magnetic equivalent circuit and a dynamical qd0 reference frame model. The needed equations for the proper design of the LSPMS motors were presented. The operating point of the permanent magnet was determined analytically, considering the armature reaction effect.

Electrical motors consist of many parts, which are mainly electrical and mechanical, such as stator and rotor winding, rotor bars or damper windings, bearings, and rotor magnets in the Permanent Magnet Synchronous Motor (PMSM) and LSPMS motors. Those parts are exposed to failure due to many reasons, such as mistakes during repair and non-ideal operational environment, such as dust, vibrations, high temperature, etc. [9,10]. The electrical faults include windings short-circuit or inter-turn faults, and rotor permanent magnets demagnetization due to armature reaction effect. Mechanical faults consist of airgap eccentricity, broken rotor bars, and bearing damage. Among those faults, 40% of faults are related to bearing damage, 38% are stator related faults, 10% are rotor related faults, and 12% are external and other faults [11]. As any other motor, LSPMS motors may experience different kinds of faults. Airgap Eccentricity related faults are responsible for about 5 to 15 % of the total faults. It happens when there is a non-uniform distance between stator inner race and rotor outer race, which in turn, will result in an unbalance magnetic field in the airgap. The unbalanced magnetic flux will lead to having a non-uniform flux linking both the stator and rotor circuits. This will affect mainly the motor inductances due to the non-uniform air gap distance. Moreover, the synchronization process and time variation of the stator current will be affected. There will be an injection of harmonics in the stator line current [12]. Eccentricity faults are attributed to many reasons, such as overloading, overstress in thermal and pressure circumstances, rotor misalignments, and bearing faults, which will result in motor shaft vibration and acoustic noise, which could completely damage the motor [13]. Three types of eccentricity can occur in electrical motors: static, dynamic, and mixed eccentricity [9]. The static eccentricity fault condition occurs when the rotor is at a constant offset from the stator symmetry axis and rotates with the new center of rotation, which has a specific offset from stator center. Thus, the minimum air gap distance is constant. However, the motor rotor is rotating at another center, rather than the symmetry center [14].

Various methods were presented in the literature for modeling and simulating eccentricity faults, the majority are based on the finite element analysis method or coupled magnetic circuit approach and the modified winding function method (MWFM). The finite element method (FEM) is considered as one of the most accurate methods to simulate the motor behavior [15,16,17,18]. However, it requires a knowledge of detailed motor parameters in addition to the internal construction of the motor to represent the geometry model. This method is a time-consuming process, even for a powerful computer, and will not provide as simple analytical modeling expressions as the MWFM method [19]. On the other hand, the MWFM can be used to calculate the motor inductances under the static eccentricity fault to represent the non-uniform airgap [20]. Moreover, the results obtained from the MWFM are verified using the finite element method for validation purposes.

The MWFM depends on the winding configuration and the positional relationship of the stator and rotor coils along with the physical motor dimensions. It is based on the Ampere and Gauss laws that drive the various motor inductances, which include the effect of space harmonics [21]. The MWFM is widely introduced to simulate and model electrical motors both in healthy and under internal motor fault conditions in both synchronous and induction machines. It is utilized for the first time in [22] for internal winding fault analysis, and consequently in [23,24], for rotor faults. LSPMS motors behave like induction motors in the starting process and operate in a steady state as the synchronous motors. Therefore, it is important to review the methods of eccentricity modeling, both in induction and synchronous motors. Nandi et al. [25] modeled the mixed eccentricity in three phase induction motors based on the MWFM along with the finite element analysis to compute the motor inductances. Frequency analysis of stator current is utilized as a detection scheme. The same methodology was followed in [26] by performing a dynamical simulation to the induction motor under the mixed eccentricity fault. In [13], the authors present a closed form expression for calculating the inductance matrices in induction motors under eccentricity conditions where a new technique based on the winding function method is proposed to model the eccentricity. Toliyat et al. [12] introduced the MWFM to present a new model of the induction motor under airgap eccentricity. The authors developed the mathematical model based on the magnetically coupled circuit approach, and represent the eccentricity by modifying the inductance matrices of the motor. Wu et al. [27] used the same principle in separating the load oscillation and rotor fault effects in the stator current based on MWFM.

Eccentricity analysis in synchronous motors is similar to the induction motors. The authors in [28] present a new index for detecting the dynamic eccentricity, which is based on the modified winding function method in the round rotor synchronous motors. This index is based on processing the developed torque using the time series data mining (TSDM) method. The authors also investigated the load variation in [29] under eccentricity conditions. The authors showed that the load variation has no impact on the amplitude of stator current harmonics since the motor is running at a constant speed. The same approach was used in [30] to present a novel method for modeling dynamic eccentricity in salient-pole synchronous motors. The coupled magnetic circuit approach has been used for simulating this fault. The main key for this simulation is to determine the machine inductances under healthy and faulty cases. Besides, the winding function method was used in [31] for calculating the mutual inductance between the stator and rotor in the LSPMS motor. The authors presented the rotor as portions; each part has a high permeability and is separated by one layer of low permeable material. The finite element method is used to confirm the simulation results.

In regards to the analysis of LSPMS motors under static eccentricity fault condition, few attempts were recently published [14,32,33]. They all simulate the motor using the FEM. Therefore, this paper is devoted for the development of a mathematical model and analyzing the performance of LSPMS motor under static eccentricity fault condition. The developed model will be used to simulate the motor transient and steady state performance using MATLAB^® software in qd0 synchronously reference frame. Results will be compared to FEM analysis using JMAG^® to judge the model accuracy. The effect of static eccentricity on the motor’s line current, speed, and torque characteristics will be investigated. The paper is organized as follows. The next section present the development of mathematical models for the LSPMS motor under healthy and static eccentricity conditions. Section 2 will also include the developed JMAG FEM model. Section 3 presents detailed simulation results and analysis for both the developed mathematical model and JMAG FEM simulation. Finally, some conclusions and possible future work is presented in Section 4.

2. Mathematical and FEM Modeling of LSPMS Motors

This section presents the dynamical and FEM modeling of LSPMS motor under healthy and static eccentricity condition.

2.1. Dynamical Model of LSPMS Motor under Healthy Conditions

The healthy abc-LSPMS motor model is well known [2,34]. To drive this model, some assumptions should be made, which are: the stator winding distribution is sinusoidal, the iron magnetic permeability is infinite, and the magnetic operating point is assumed to have a fixed value and located in the linear region. Figure 1 shows the motor stator and rotor circuits assuming two pole motor. The origin of coordinate $O$, and the winding axis of stator and rotor, as well as the quadrature and direct axis of magnet are also shown on the same figure.

By using the coupled magnetic circuit approach and Kirchhoff’s law, the stator and rotor voltage equations can be written as:

$$\hspace{1em}\left[\begin{array}{c}{V}_{s,abc}\\ V_{r,abc}\end{array}\right]=\left[\begin{array}{cc}{r}_{s,abc}& 0\\ 0& r_{r,abc}\end{array}\right]\left[\begin{array}{c}{i}_{s,abc}\\ i_{r,abc}\end{array}\right]+\frac{d}{dt}\left[\begin{array}{c}{\lambda }_{s,abc}\\ {\lambda }_{r,abc}\end{array}\right]$$

(1)

where $V_{s,abc}$ and $V_{r,abc}$ are the stator and rotor terminal voltages. $r_{s,abc}$ and $r_{r,abc}$ are the stator and rotor phase resistances. $i_{s,abc}$ and $i_{r,abs}$ are stator and rotor phase currents. ${\lambda }_{s,abc}$ and ${\lambda }_{r,abc}$ are the stator and rotor flux linkages. For conciseness, matrix expressions are used to donate three phase variables, as follows.

$$\begin{gathered} V_{s,abc}={\left[ V_{sa}\hspace{1em}{V}_{sb}\hspace{1em}{V}_{sc} \right]}^T,\hspace{1em}{i}_{s,abc}={\left[ i_{sa}\hspace{1em}{i}_{sb}\hspace{1em}{i}_{sc} \right]}^T \\ {V}_{r,abc}={\left[ V_{ra}\hspace{1em}{V}_{rb}\hspace{1em}{V}_{rc} \right]}^T,\hspace{1em}{i}_{r,abc}={\left[ i_{ra}\hspace{1em}{i}_{rb}\hspace{1em}{i}_{rc} \right]}^T \\ {\lambda }_{s,abc}={\left[ {\lambda }_{sa}\hspace{1em}{\lambda }_{sb}\hspace{1em}{\lambda }_{sc} \right]}^T,\hspace{1em}{\lambda }_{r,abc}={\left[ {\lambda }_{ra}\hspace{1em}{\lambda }_{rb}\hspace{1em}{\lambda }_{rc} \right]}^T \\ {r}_{s,abc}=\left[\begin{array}{ccc}{r}_s& 0& 0\\ 0& r_s& 0\\ 0& 0& r_s\end{array}\right],\hspace{1em} r_{r,abc}=\left[\begin{array}{ccc}{r}_r& 0& 0\\ 0& r_r& 0\\ 0& 0& r_r\end{array}\right] \end{gathered}$$

(2)

The stator magnetic fluxes ${\lambda }_{s,abc}$ are generated by three different sources, which are: the flux created by self-inductances $L_{ss}$, the second is formed by the mutual inductance between stator and rotor $L_{sr}$, and finally, the permanent magnet flux ${\lambda }_{s,pm}$. The same analogy is applied for the rotor magnetic flux ${\lambda }_{r,abc}$. The total flux linkage equations can be written in matrix format as:

$$\left[\begin{array}{c}{\lambda }_{s,abc}\\ {\lambda }_{r,abc}\end{array}\right]=\left[\begin{array}{cc}{L}_{ss}& L_{sr}\\ L_{rs}& L_{rr}\end{array}\right]\left[\begin{array}{c}{i}_{s,abc}\\ i_{r,abc}\end{array}\right]+\left[\begin{array}{c}{\lambda }_{s,pm}\\ {\lambda }_{r,pm}\end{array}\right]$$

(3)

$$\begin{gathered} L_{ss}=\left[\begin{array}{ccc}{L}_{asas}& L_{asbs}& L_{ascs}\\ L_{bsas}& L_{bsbs}& L_{bscs}\\ L_{csas}& L_{csbs}& L_{cscs}\end{array}\right], L_{rr}=\left[\begin{array}{ccc}{L}_{arar}& L_{arbr}& L_{arcr}\\ L_{brar}& L_{brbr}& L_{brcr}\\ L_{crar}& L_{crbr}& L_{crcr}\end{array}\right] \\ {L}_{sr}=\left[\begin{array}{ccc}{L}_{asar}& L_{asbr}& L_{ascr}\\ L_{bsar}& L_{bsbr}& L_{bscr}\\ L_{csar}& L_{csbr}& L_{cscr}\end{array}\right] , L_{rs}={\left(L_{sr}\right)}^T\hspace{1em} \end{gathered}$$

(4)

$$\hspace{1em}{\lambda }_{s,pm}={\lambda }_m\left[\begin{array}{c}\sin \left({\theta }_r\right)\\ \sin \left({\theta }_r-\frac{2\pi }{3}\right)\\ \sin \left({\theta }_r+\frac{2\pi }{3}\right)\end{array}\right] , {\lambda }_{r,pm}={\lambda }_m\left[\begin{array}{c}0\\ \frac{\sqrt{3}}{2}\\ -\frac{\sqrt{3}}{2}\end{array}\right]\hspace{1em}$$

(5)

where $L_{ss},\hspace{1em}{L}_{sr},\hspace{1em}{L}_{rs},\hspace{1em}{L}_{rr}$ are the stator self-inductance, stator to rotor mutual inductance, rotor to stator mutual inductance, and rotor self-inductance matrices, respectively. ${\lambda }_{s,pm}$ and ${\lambda }_{r,pm}$ are the permanent magnet flux component for stator and rotor, respectively. ${\lambda }_m$ is the peak value of the PM flux. ${\theta }_r$ is the reference angle of the PM flux vector. The diagonal elements of $L_{ss}$ and $L_{rr}$ are the self-inductances of stator and rotor windings, while the off-diagonal elements are the mutual inductances. They can be calculated under healthy condition using the winding function method approach [28]. Park transformation [35] in the rotor reference frame is used to transform the abc-LSPMS motor model into qd0-model. This yields (6)–(10), which complete the model for the healthy LSPMS motor.

$$\left[\begin{array}{c}{\lambda }_{qs}\\ {\lambda }_{ds}\\ {\lambda }_{0s}\\ {\lambda }_{qr}^{\prime }\\ {\lambda }_{dr}^{\prime }\\ {\lambda }_{0r}^{\prime }\end{array}\right]=\left[\begin{array}{cccccc}{L}_{ls}+L_{mq}& 0& 0& L_{mq}& 0& 0\\ 0& L_{ls}+L_{md}& 0& 0& L_{md}& 0\\ 0& 0& L_{ls}& 0& 0& 0\\ L_{mq}& 0& 0& L_{lr}+L_{mq}& 0& 0\\ 0& L_{md}& 0& 0& L_{lr}+L_{md}& 0\\ 0& 0& 0& 0& 0& L_{lr}\end{array}\right]\left[\begin{array}{c}{i}_{sq}\\ i_{sd}\\ \begin{array}{c}{i}_{s0}\\ \begin{array}{c}{i}_{rq}^{\prime }\\ i_{rd}^{\prime }\\ i_{r0}^{\prime }\end{array}\end{array}\end{array}\right]+\left[\begin{array}{c}0\\ {\lambda }_m\\ 0\\ 0\\ {\lambda }_m^{\prime }\\ 0\end{array}\right]$$

(6)

$$V_{s,qd0}=r_s{i}_{s,qd0}+\omega \left[\begin{array}{ccc}0& 1& 0\\ -1& 0& 0\\ 0& 0& 0\end{array}\right]{\lambda }_{s,qd0}+\frac{d}{dt}{\lambda }_{s,qd0}$$

(7)

$$V_{r,qd0}=r_r^{\prime }{i}_{r,qd0}^{\prime }+\frac{d}{dt}{\lambda }_{r,qd0}^{\prime }$$

(8)

where $\omega$ is the electrical angular speed. $L_{mq}$ and $L_{md}$ are the qd0 magnetizing inductances, respectively. $L_{ls},L_{lr}$ are the leakage inductances of stator and rotor windings, respectively. ${\lambda }_m^{\prime }=i_m^{\prime }{L}_{md}$ and $i_m^{\prime }$ is the equivalent magnetization current referred to stator side. The prime notations in the last equations indicates a variable being referred to stator side. It should be noticed that the rotating reference frame is assumed to be on the rotor side. Therefore, the permanent magnet flux component, at the stator side, is changed as a sinusoidal wave form. In qd0 reference frame, this component is available along the d-axis only for the stator and rotor.

The mechanical speed is expressed in (9) as:

$$\hspace{1em}{\omega }_r\left(t\right)=\frac{P}{2J}{{\displaystyle \int }}^{\text{}}\left(T_{em}+T_{mech}-T_{damp}\right)dt$$

(9)

where $P$ is the number of poles, $J$ is the motor inerita, $T_{mech}$ is the input mechanical load, $T_{damp}$ is the damping torque, and ${\omega }_r$ is the rotor angular velocity. Finally, the electromechanically torque that is developed by the motor can be obtained from the input power being transferred into the airgap [34], which is given by:

$$T_{em}=\frac{3}{2}\frac{P}{2}\left( {\lambda }_{ds}{i}_{qs}-{\lambda }_{qs}{i}_{ds} \right)$$

(10)

2.2. Dynamical Model of LSPMS Motor under Static Eccentricity Condition

Eccentricity, in general, will affect the motor inductances. Thus, the previously presented inductance matrices in a healthy condition, as shown in (4), are not valid for the model under the static eccentricity fault condition [36]. In this paper, coupled magnetic circuit approach is used to represent the abc-LSPMS motor model under static eccentricity fault condition, while the MWFM is used for estimating the motor inductances. This section will include the modeling of the airgap function and inductance calculation under this fault.

2.2.1. Modeling of Static Eccentricity

Eccentricity is a mechanical related fault; the occurrence of static eccentricity will change the motor rotational axis symmetry. Therefore, considering the non-uniform airgap distribution is necessary. According to [36], the airgap permeance is inversely proportional to the airgap length, as shown in (11). Thus, the airgap permeance can be neglected for points far from the rotor poles.

$$P\left(\varphi \right)={\mu }_0\frac{r_{av}\left(\varphi \right)}{g\left(\varphi \right)}\hspace{1em}$$

(11)

where ${\mu }_0,\hspace{1em}{r}_{av}\left(\varphi \right),$ $g\left(\varphi \right)$, and $\varphi$ are the permeability in space, the average radius of airgap, the airgap function, and the base angle for measurement, respectively. Figure 2 presents the motor geometry with static eccentricity. The rotor center of rotation $o^{\prime }$ is offset from the center of stator $o$ to represent an eccentric rotor. Vector $\overrightarrow{ O{O}^{\prime }}$ is called the static eccentricity vector. The degree of static eccentricity ${\delta }_s$ is defined as the ratio between the magnitude of static eccentricity vector $\overrightarrow{ O{O}^{\prime }}$ to the symmetrical airgap length $g_0$ [14], yields:

$${\delta }_s=\frac{\left|\overrightarrow{ O{O}^{\prime }}\right|}{g_0}$$

(12)

This vector is fixed in space under static eccentricity and its position is characterized by angle ${\theta }_0$. This angle describes the stationary position of the minimum airgap along the inner stator surface. For an arbitrary point M located on the rotor outer race, the distance between this point and the rotor center is equal to rotor radius [37], or:

$$OM={\delta }_{s}g\cos \left(\varphi -{\theta }_0\right)+\sqrt{R_r^2-{\delta }_s^2{g}_0^2{\sin }^2\left(\varphi -{\theta }_0\right)}\hspace{1em}$$

(13)

where $R_r$ is the rotor outer radius. For stator inner radius of $R_s$, the airgap function can be obtained as:

$$g\left(\varphi \right)=R_s-OM\hspace{1em}$$

(14)

Hence, the rotor radius is much larger than the airgap length, the second term in (13) can be represented by $R_r$. Therefore, (14) can be written as:

$$g\left(\varphi \right)=g_0\left[1-{\delta }_s\cos \left(\varphi -{\theta }_0\right)\right]$$

(15)

For the calculation of the inductances in the MWFM, the inverse of the airgap function is needed. It is represented in (16) as [37]:

$$g^{-1}\left(\varphi \right)=\frac{1}{g_0\left[1-{\delta }_s\cos \left(\varphi -{\theta }_0\right)\right]}\hspace{1em}$$

(16)

The above equation can be expressed in Fourier series form as:

$$g^{-1}\left(\varphi \right)=A_0+{\displaystyle \sum }_n^{\infty }{A}_n\cos n\left(\varphi -{\theta }_0\right)\hspace{1em}n=1,2..\hspace{1em}$$

(17)

The series coefficients $A_0$ and $A_n$ can be obtained as:

$$A_0=\frac{1}{g_0\sqrt{1-{\delta }_s^2}} ;\hspace{1em}{A}_n=\frac{2}{g_0\sqrt{1-{\delta }_s^2}}{\left[\frac{1-\sqrt{1-{\delta }_s^2}}{{\delta }_s^2}\right]}^n\hspace{1em}$$

(18)

The inverse function for the first term can be written as:

$$g^{-1}\left(\varphi \right)=\frac{1}{g_0\sqrt{1-{\delta }_s^2}}+2\frac{1-\sqrt{1-{\delta }_s^2}}{g_0{\delta }_s^2\sqrt{1-{\delta }_s^2}}\cos \left(\varphi -{\theta }_0\right)$$

(19)

2.2.2. Inductance Calculations under Static Eccentricity

The MWFM will be used to derive the motor inductances under static eccentricity fault. A detailed discussion about this method is reported in [20]. This method assumes two set of arbitrary winding A and B. The self and mutual inductance between these two winding can be written as:

$$L_{aa}={\mu }_{0}rl {{\displaystyle \int }}_0^{2\pi }{N}_a\left(\varphi \right) n_a\left(\varphi \right)g^{-1}\left(\varphi \right) d\varphi \hspace{1em}$$

(20)

$$L_{ba}={\mu }_{0}rl {{\displaystyle \int }}_0^{2\pi }{N}_a\left(\varphi \right) n_b\left(\varphi \right)g^{-1}\left(\varphi \right) d\varphi \hspace{1em}$$

(21)

where $L_{aa}$ is the self-inductance of winding a, $L_{ba}$ is the mutual inductance linking winding b due to current flow in winding a. $\varphi$ is a reference angle along the stator inner surface. The terms $N_a\left(\varphi \right)$, $n_a\left(\varphi \right)$, and $n_b\left(\varphi \right)$ are the modified winding function of winding a, and the turn functions for winding a and b, respectively. The modified winding function for arbitrary winding a is defined as:

$$\hspace{1em}{N}_a\left(\varphi \right)=n_a\left(\varphi \right)-{{\displaystyle \int }}_0^{2\pi }\frac{n_a\left(\varphi \right)g^{-1}\left(\varphi \right)}{2\pi g_{\mathrm{av}}^{-1}\left(\varphi \right)}d\varphi$$

(22)

where $g_{av}^{-1}\left(\varphi \right)$ is the average value of (19) and equal to:

$$g_{av}^{-1}\left(\varphi \right)=A_0$$

(23)

The inductance matrices as presented in (4), for the healthy condition, should be recalculated under static eccentricity fault condition. For this purpose, MATLAB^® symbolic integration method is used to calculate the total 27 inductance elements in the matrices based on (20) and (21). The flux linkage equations under static eccentricity fault condition can be written as:

$$\left[\begin{array}{c}{\lambda }_{s,abc}\\ {\lambda }_{r,abc}\end{array}\right]=\left[\begin{array}{cc}{L}_{ss}^{SE}& L_{sr}^{SE}\\ L_{rs}^{SE}& L_{rr}^{SE}\end{array}\right]\left[\begin{array}{c}{i}_{s,abc}\\ i_{r,abc}\end{array}\right]+\left[\begin{array}{c}{\lambda }_{s,pm}\\ {\lambda }_{r,pm}\end{array}\right]\hspace{1em}$$

(24)

$$\begin{gathered} L_{ss}^{SE}=\left[\begin{array}{ccc}{L}_{asas}^{SE}& L_{asbs}^{SE}& L_{ascs}^{SE}\\ L_{bsas}^{SE}& L_{bsbs}^{SE}& L_{bscs}^{SE}\\ L_{csas}^{SE}& L_{csbs}^{SE}& L_{cscs}^{SE}\end{array}\right] ,L_{rr}^{SE}=\left[\begin{array}{ccc}{L}_{arar}^{SE}& L_{arbr}^{SE}& L_{arcr}^{SE}\\ L_{brar}^{SE}& L_{brbr}^{SE}& L_{brcr}^{SE}\\ L_{crar}^{SE}& L_{crbr}^{SE}& L_{crcr}^{SE}\end{array}\right] \\ {L}_{sr}^{SE}=\left[\begin{array}{ccc}{L}_{asar}^{SE}& L_{asbr}^{SE}& L_{ascr}^{SE}\\ L_{bsar}^{SE}& L_{bsbr}^{SE}& L_{bscr}^{SE}\\ L_{csar}^{SE}& L_{csbr}^{SE}& L_{cscr}^{SE}\end{array}\right], L_{rs}^{SE}={\left(L_{sr}^{SE}\right)}^T \end{gathered}$$

(25)

where $L_{ss}^{SE}$, $L_{sr}^{SE}$, $L_{rs}^{SE}$, $L_{rr}^{SE}$ are the inductance matrices under the static eccentricity fault condition. The diagonal elements of $L_{ss}^{SE}$ and $L_{rr}^{SE}$ in (25) are the self-inductances of stator and rotor windings, and they are calculated using (20), while the off-diagonal elements are the mutual inductances and they are calculated using (21). Park transformation [35] in the rotor reference frame is used to transform the abc-LSPMS motor flux linkage equations in (24) into qd0-reference frame. This yield Equations (26)–(30), which complete the model for the faulty LSPMS motor under static eccentricity.

$$\left[\begin{array}{c}{\lambda }_{sq}\\ {\lambda }_{sd}\\ {\lambda }_{s0}\\ {\lambda }_{rq}^{\prime }\\ {\lambda }_{rd}^{\prime }\\ {\lambda }_{r0}^{\prime }\end{array}\right]=\left[\begin{array}{cc}{L}_{s,qd0}^{SE}& L_{sr,qd0}^{SE}\\ L^{SE}{}_{rs,qd0}& L_{r,qd0}^{SE}\end{array}\right]\left[\begin{array}{c}{i}_{sq}\\ i_{sd}\\ \begin{array}{c}{i}_{s0}\\ \begin{array}{c}{i}_{rq}^{\prime }\\ i_{rd}^{\prime }\\ i_{r0}^{\prime }\end{array}\end{array}\end{array}\right]+\left[\begin{array}{c}0\\ {\lambda }_m\\ 0\\ 0\\ {\lambda }_m^{\prime }\\ 0\end{array}\right]\hspace{1em}$$

(26)

$$L_{s,qd0}^{SE}=\left[\begin{array}{ccc} \frac{3}{2}(L_{ss,q}+L_e{\cos }^2\left({\theta }_r-{\theta }_0\right))+L_{ls}& \frac{3}{4}{L}_e\sin 2\left({\mathsf{\theta }}_{\mathrm{r}}-{\theta }_0\right)& 0\\ \frac{3}{4}{L}_e\sin 2\left({\theta }_r-{\theta }_0\right)& \frac{3}{2}(L_{ss,d}+L_e{\sin }^2\left({\theta }_r-{\theta }_0\right))+L_{ls}& 0\\ 0& 0& L_{ls}\end{array}\right]$$

(27)

$$L_{sr,qd0}^{SE}=\left[\begin{array}{ccc}\frac{3}{2}(L_{sr,q}+L_e{\cos }^2\left({\theta }_r-{\theta }_0\right))& \frac{3}{4}{L}_e\sin 2\left({\theta }_r-{\theta }_0\right)& 0\\ \frac{3}{4}{L}_e\sin 2\left({\theta }_r-{\theta }_0\right)& \frac{3}{2}(L_{sr,d}+L_e{\sin }^2\left({\theta }_r-{\theta }_0\right))& 0\\ 0& 0& 0\end{array}\right]$$

(28)

$$\hspace{1em}{L}^{SE}{}_{rs,qd0}=L_{sr,qd0}^T$$

(29)

$$L_{r,qd0}^{SE}=\left[\begin{array}{ccc}\frac{3}{2}(L_{rr,q}+L_e{\cos }^2\left({\theta }_r-{\theta }_0\right))+L_{lr}& \frac{3}{4}{L}_e\sin 2\left({\theta }_r-{\theta }_0\right)& 0\\ \frac{3}{4}{L}_e\sin 2\left({\theta }_r-{\theta }_0\right)& \frac{3}{2}(L_{rr,d}+L_e{\sin }^2\left({\theta }_r-{\theta }_0\right))+L_{lr}& 0\\ 0& 0& L_{lr}\end{array}\right]$$

(30)

$$L_e=- \frac{\pi {\mu }_{0}lr}{8g_0}\frac{A_1^2}{A_0} {\left[\frac{N}{P_p}\right]}^2$$

(31)

$$L_{ss,q}=L_{rr,q}=\frac{\pi {\mu }_{0}lr}{4g_0}{\left[\frac{N}{P_p}\right]}^2, L_{ss,d}=L_{rr,d}=\frac{\pi {\mu }_{0}lr}{4g_0\zeta }{\left[\frac{N}{P_p}\right]}^2$$

(32)

where $L_{s,qd0}^{SE}$, $L_{sr,qd0}^{SE}$, $L_{rs,qd0}^{SE}$, and $L_{r,qd0}^{SE}$ are the transformed stator self-inductance, stator to rotor mutual inductance, rotor to stator mutual inductance, rotor self-inductance matrices under the faulty condition, respectively. $L_e$ is an inductance component due to static eccentricity. $r$ is the average radius of the airgap. $l$ is the average length of the motor stack. $P_p$ is the number of pole pairs. $\zeta$ is the saliency ratio. $N$ is the equivalent number of winding turns.

In conclusion, it should be noticed that static eccentricity have no effects on the form of the voltage equations obtained in (7) and (8). Therefore, the LSPMS motor mathematical model under eccentricity condition is represented by (7)–(9) and (26)–(32).

2.3. Modeling of LSPMS Motor Using FEM

The finite element method is a numerically based approach used for computing the motor performance characteristics and parameters, such as current, torque, speed, flux density, and inductances. This method will handle the complex motor geometry, consider the mixed set of materials that are involved in the motor parts, and, in addition, monitor the non-linear behavior of these materials such as permanent magnets. In this method, the motor geometry will be divided into an enormous number of elements joined with nodes by using a proper discretization and mesh generation algorithms. A set of equations represents each element, and the entire set of equations and their solutions will implement the motor model.

To ensure the accuracy of developed mathematical model, FEM based analysis is performed for three phase LSPMS motor (1-hp, 415 V, Y-connected, four-pole) utilizing JMAG^® designer software. JMAG is an electromagnetic simulation program that was developed by JSOL Corporation. JMAG provides a two-dimensional (2-D)/three-dimensional (3-D) platform to construct motor models. Material properties can be assigned to each part of the model. JMAG includes a material library based on commercial and industrial materials. The material properties and names are obtained from the manufacturer’s data sheets and brochures. Additional materials can be easily added, along with their complete mechanical, electrical and magnetic properties. Moreover, the available material can be customized based on the application. Figure 3 shows the motor that was used in the study. To create the finite element model, the motor geometry was drawn using Auto-CAD^TM and then imported into JMAG^TM.

The motor geometry specifications and electrical parameters are shown in

Table 1. LSPMS Motor Specifications.

Parameter	Value
Rated power (W)	750
Rated Voltage (V)	415
Stator phase resistance (Ohm)	19.15
Stator leakage inductance (mH)	0.004
Number of poles	4
Frequency (Hz)	50
Air-gap length (mm)	0.3
Outer/inner stator diameter (mm)	120/75
Number of stator/rotor slots	24/16
Axial length of stator core (mm)	75
Number of turns per slot	139
Height of stator yoke (mm)	45
Height of stator/rotor slots (mm)	13/9.5
Magnet material	Recoma-24HE
Remanent of magnet (T)	1.02

. Figure 4 presents the geometry model of the motor. The model consists of all the geometrical parts, which include the laminated silicon steel core of the stator and rotor, stator copper windings, shaft, rotor casts aluminum bars, and the permanent magnet of type 24HE. The monolayer topology is used for stator winding distribution, and this will reduce the total harmonics in the line voltage and current. The transient iterative solver of JMAG^®, based on incomplete Cholesky conjugate gradient (ICCG) method, is used to solve for magnetic field distribution in the model, as well as the other motor performance characteristics, such as stator current, speed, and torque.

The 2-D mesh model of the motor is created with 70,729 nodes and 133,803 elements. Figure 5 shows the resulting mesh model of LSPMS motor.

3. Simulation Results

The simulation results are divided into two main parts. First, a comparison between the simulation results of the MATLAB^® and JMAG^® models is presented under healthy conditions. The second part will contain the simulation results under static eccentricity fault conditions. This section will include a simulation and comparison of 21 cases, using JMAG^® and MATLAB^®. They will cover a wide range of eccentricity degrees under different loading conditions. Also, the angle effect of the static eccentricity vector along the motor circumference has been studied. The simulation results will be compared in terms of stator current, speed, and torque characteristics. Figure 6 presents the flow chart of modeling and simulation using MATLAB^®.

The simulation starts by determining the inductance matrices given in (26). Simulation is initialized for the 1-hp LSPMS motor with parameters given in Table 1. Three phase sinusoidal voltages are applied to the motor terminals. By rearranging (7) and (8), the quadrature and direct fluxes are calculated after performing the integration. The inductance matrices calculated in (26) will be used to determine the current in qd0 reference frame. From the fluxes and currents calculated, the electromagnetic torque, rotor angular position ${\theta }_r$, and speed ${\omega }_r$ are determined using (9). The performance of the LSPMS motor will be investigated under the normal operating loads. The selected degrees of static eccentricity (SE) based on (12) will be in the range of 13 to 33%. Under 10% SE, a small variation is expected in the motor line current, speed, and torque characteristics. This is due to intrinsic SE in the normal operation of the motor, which was estimated to be around the 10% [32].

To verify the accuracy of the derived model under different loading and eccentricity conditions, the LSPMS motor is simulated using JMAG^® and MATLAB^®. Figure 7, Figure 8 and Figure 9 show the stator current, electromagnetic torque, and speed for a load of 1.5 N·m. Results show a good agreement between JMAG^® and MATLAB^® results from the developed model. The stator current shows very similar time variations in transient and steady state operation. Moreover, the torque and speed characteristics of the developed mathematical model follow the same behavior of the FEM simulation results. It should be noticed that the inductance matrix is calculated using the MWFM. The inductance matrix is not constant over the entire simulation time. In transient region, the inductance matrix is changing due to current as well as flux changes. The simulation is started by setting some initial values for the flux. Therefore, to find the corresponded current value due to that flux, the inductance matrix is derived as a function of the rotor angular position. Thus, the calculated inductance value will be a function of that current, torque, and angular position. In steady state, the current and torque quantities will reach a steady state value. Therefore, the inductance matrix elements will have a constant value. This applies for both the healthy and faulty conditions.

To demonstrate the accuracy of the developed model under eccentricity condition, Figure 10 depicts the stator current under 19% static eccentricity for 2.1 N·m load. The similarity of time variations between MATLAB^® and FEM results are clear.

Figure 11 shows a comparison between the healthy and faulty stator currents under 19% SE and a load of 2.1 N·m. It is quite clear that under this degree of SE, it takes a longer time to attain a steady state value than a healthy condition Also, higher current peaks are presented in the transient region.

For further comparison between FEM and the developed mathematical model results,

Table 2. The percentage correlation and normalized mean square error (NMSE) of simulated cases.

Case Number	Load N·m	%${\mathit{\delta }}_{\mathit{s}}$	%Correlation	%NMSE
Case 1	0	0	82.3	0.8
Case 2		16	83.1	0.7
Case 3	1.5	0	81.5	2.1
Case 4		16	83.7	2.0
Case 5		22	80.3	2.3
Case 6	2	0	81.6	2.4
Case 7		16	80.0	2.4
Case 8	2.1	0	76.0	2.6
Case 9		13	76.3	3.1
Case 10		19	75.0	3.5
Case 11		30	76.2	3.2
Case 12	2.2	0	77.4	3.2
Case 13		13	76.8	3.1
Case 14		16	76.0	3.8
Case 15	2.3	0	80.1	3.2
Case 16		18	74.1	4.1
Case 17		30	79.0	3.4
Case 18		33	77.1	3.7
Case 19	2.4	0	75.0	3.5
Case 20		16	74.0	4.1
Case 21		33	76.0	3.2

presents the correlation and the normalized mean square error (NMSE) for 21 cases of different loading conditions and static eccentricity degrees. As reported in the literature, it is useful to pay more attention to the stator current as the fault diagnostic techniques, such as motor current signature analysis (MCSA), are based on the stator current waveform. The NMSE is defined as:

$$\%NMSE=\frac{{{\displaystyle \sum }}_{i=1}^N{\left( \left(\frac{X_i}{Max\left(x\right)}\right)- \left(\frac{X_i^{\prime }}{Max\left(x^{\prime }\right)}\right)\right)}^2}{N} \times 100 \%$$

(33)

where $X_i/Max\left(x\right)$ and $X_i^{\prime }/Max\left(x^{\prime }\right)$ are the JMAG^® and MATLAB^® stator current time variation values divided by the maximum value of each time series, respectively. $N$ is the number of samples per each series.

The correlation coefficient is a statistical measure of dependence between two variables, and it is also an indicator of how two variables fluctuate together.

Under different levels of static eccentricity (0% SE, 13% SE, 16% SE, and 19% SE) and a load of 2.1 N·m, the stator current and speed characteristics have different behaviors, Figure 12 and Figure 13. It is clear that as the eccentricity level increases, the synchronization time is also increased. A slight difference is shown in the results between the healthy and 13% static eccentricity condition. However, more peaks, ripples, and longer synchronization times are observed for 16% and 19% static eccentricity before reaching the steady state region.

It is worth mentioning that the use of this accurate mathematical model for developing an online motor current signature diagnostic system can be integrated to the electric utilities and industrial plants predictive maintenance systems in order to run their facilities at high efficiency and to reduce maintenance and failure costs.

Sensitivity analysis is quite important to highlight the main sources of uncertainty and disturbances, which affects the simulation results. As such, and in order to test the robustness and sensitivity of the developed mathematical model, some parameters have been selected and varied. For example, the stator resistance $r_s$ and rotor equivalent resistance $r_r$ values are varied by ±10% of the original parameter value. The simulation is performed for 36 cases, which represent different loading conditions and eccentricity levels. The percentage normalized mean square error (%NMSE) is calculated to show the difference in results corresponding to the original and the modified parameter values. This comparison is conducted for the motor line current characteristic, which is the main output characteristic in this paper.

Table 3. The percentage NMSE of simulated cases.

Eccentricity Percentage	%NMSE
	Under Healthy state (${\delta }_s=0$)			Under ${\delta }_s=16\%$			Under ${\delta }_s=22\%$
Load Parameter	0 N·m	1.5 N·m	2.4 N·m	0 N·m	1.5 N·m	2.4 N·m	0 N·m	1.5 N·m	2.4 N·m
$+10\% of r_r$	0.20	2.03	4.62	0.99	2.70	4.43	0.99	2.63	4.39
$-10\% of r_r$	0.65	1.47	3.28	0.24	2.25	2.45	1.24	2.49	9.40
$+10\% of r_s$	1.50	2.93	10.00	0.76	3.03	11.48	1.18	3.67	18.55
$-10\% of r_s$	2.59	2.81	10.31	1.76	3.79	11.45	1.15	3.35	13.71

presents the percentage NMSE for each simulated case. The results show clearly the robustness of the proposed mathematical model against variation of these motor parameters.

4. Conclusions

This paper proposes a new mathematical and simulation models of LSPMS motor under static eccentricity fault conditions. The simulation is performed to investigate and compare the motor performance under healthy and static eccentricity fault conditions while considering different operating loads. The results obtained for the developed mathematical model under static eccentricity conditions are in good agreement with the FEM simulation results. Besides, results show that the effect of static eccentricity on the stator current depends on the degree of eccentricity, which will affect the motor synchronization time and the time variation of stator current in transient period with less significance on the steady state performance. The results are in good agreement with the findings reported in the literature when using the finite element method. This confirms the precision of the proposed mathematical model for eccentricity conditions. On the other hand, the problem of the eccentricity fault in LSPMS motors is far from being completed by only relying on a limited number of case studies and software packages. Therefore, the analysis could be extended to include the dynamic and mixed eccentricity fault conditions. In addition, the mathematical model can be further improved by considering the nonlinearity between the motor parameters and the output characteristics.