Comparison of nonlinear solution methods for magnetic equivalent circuits of saturated induction motors

Abstract

This work compares three nonlinear solution methods for the performance of an induction motor’s magnetic equivalent circuit model with magnetic saturation. The interrelation between magnetic flux density and permeability introduces nonlinearities in the differential system of equations. Three popular nonlinear solution methods are selected for comparison, namely (i) the Gauss–Seidel method, (ii) the Newton–Raphson method and (iii) the inverse Broyden’s method. While all three methods have been applied in this context before, no comparison study has been published to the authors’ best knowledge. The study finds that the inverse Broyden’s method is most performant in terms of the number of required iterations, the computation time per iteration and the resulting total computation time. However, for substantial saturation levels, the authors recommend a hybrid implementation of multiple solution methods to obtain robust and reliable convergence.

Keywords

Magnetic equivalent circuit nonlinear induction motor saturation Gauss–Seidel Newton–Raphson inverse Broyden

1. Introduction

Magnetic Equivalent Circuit (MEC) models can generate accurate results in a computationally efficient way when compared to e.g. finite element method (FEM) models or direct-quadrature (DQ) models [1]. In practice, magnetic saturation of the stator or rotor core can easily happen as they are designed to operate at the knee point of their magnetisation curves [2]. Hence, loading the motor above the rated torque or rotor eccentricity can cause saturation. This effect is governed by the nonlinear magnetisation of the material, which decreases material permeability when the magnetic field strength is increased. Because of this, an iterative solution method is required. Most often, the Gauss–Seidel method (GSM), also called the method of successive displacement, is used due to its straightforward implementation and acceptable rate of convergence. Such is the case in the works of e.g. [3] for a MEC of an induction motor [4], for a MEC of a permanent magnet synchronous motor [5], for a quasi 3D MEC of a double-sided permanent magnet synchronous machine and [6] for a hybrid MEC of general electrical machines. Often, the convergence rate benefits from the integration of the derivatives of the governing equations, expressed in a Jacobian matrix. This is the basis for the Newton–Raphson method (NRM), also called Newton’s method, which is used in the works of e.g. [7] for a 3D MEC of an induction motor and [8] for a hardware-in-the-loop MEC of an induction motor. As MEC’s main advantage is computational speed and as magnetic saturation contributes significantly to the total computation time, more advanced nonlinear solution methods should be explored to improve the convergence rate. Broyden’s method is a modification of NRM, which reduces the need to evaluate the exact Jacobian. In the context of MECs of induction motors, this was applied in the work of [9] and mentioned by [10] and [11]. A modified version of Broyden’s method, namely Inverse Broyden’s Method (IBM) is used in this work [12]. This article provides a performance comparison of the aforementioned solution methods for MECs of induction motors with saturation induced by overloading and rotor eccentricity. The emphasis of the analysis is on computational speed, in order to retain the strength of MEC theory even in the presence of saturation.

2. Nonlinear system modelling

2.1. Magnetic equivalent circuit

The nonlinear MEC is based on [13] and includes the effects of slotting, rotor skewing and radial rotor eccentricity. As shown in Fig. 1, the MEC couples the stator and rotor electrical circuits using a magnetic circuit that spans over the stator, air gap and rotor. The magnetic circuits connect to the electrical circuits through electromotive forces (EMFs) induced by magnetic flux changes and through magnetomotive forces (MMFs) induced by currents [3]. The magnetic circuit’s nodes are interconnected by reluctances in the stator and rotor and permeances in the air gap. The coupled behaviour is described $\begin{eqnarray}\displaystyle \hspace{-10.0pt}\left[\begin{array}{@{}cccccccc@{}}\mathbf{0} & \mathbf{0} & \mathbf{0} & {\displaystyle \frac{d}{dt}}\boldsymbol{N}_{abc}^{\intercal } & \mathbf{0} & \mathbf{0} & \boldsymbol{R}^{st} & \mathbf{0}\\ \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} & \frac{d}{dt}\boldsymbol{S}\boldsymbol{u}\boldsymbol{b}\boldsymbol{B}_{N_{rot}} & \mathbf{0} & \mathbf{0} & -\boldsymbol{R}^{rot}\\ \boldsymbol{S}\boldsymbol{u}\boldsymbol{b}\boldsymbol{F}_{N_{st}} & \mathbf{0} & \boldsymbol{{\mathcal{R}}}_{t}^{st} & -\boldsymbol{{\mathcal{R}}}_{y}^{st} & \mathbf{0} & \mathbf{0} & \boldsymbol{N}_{abc} & \mathbf{0}\\ \mathbf{0} & \mathbf{0} & -\boldsymbol{I}_{N_{st}} & \boldsymbol{S}\boldsymbol{u}\boldsymbol{b}\boldsymbol{B}_{N_{st}} & \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0}\\ \mathbf{0} & \boldsymbol{S}\boldsymbol{u}\boldsymbol{b}\boldsymbol{F}_{N_{st}} & \mathbf{0} & \mathbf{0} & \boldsymbol{{\mathcal{R}}}_{t}^{rot} & -\boldsymbol{{\mathcal{R}}}_{y}^{rot} & \mathbf{0} & \boldsymbol{I}_{N_{rot}}\\ \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} & \boldsymbol{I}_{N_{rot}} & \boldsymbol{S}\boldsymbol{u}\boldsymbol{b}\boldsymbol{B}_{N_{rot}} & \mathbf{0} & \mathbf{0}\\ \boldsymbol{{\mathcal{P}}}_{rot}^{st} & -\boldsymbol{{\mathcal{P}}}^{air} & \boldsymbol{I}_{N_{st}} & \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0}\\ -\boldsymbol{{\mathcal{P}}}^{air\intercal } & \boldsymbol{{\mathcal{P}}}_{st}^{rot} & \mathbf{0} & \mathbf{0} & -\boldsymbol{I}_{N_{rot}} & \mathbf{0} & \mathbf{0} & \mathbf{0}\\ \boldsymbol{s}\boldsymbol{e}\boldsymbol{t}_{1,N_{st}} & \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0}\\ \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{1} & \mathbf{0}\end{array}\right]\left[\begin{array}{@{}c@{}}\boldsymbol{{\psi}}^{st}\\ \boldsymbol{{\psi}}^{rot}\\ \boldsymbol{{\phi}}_{t}^{st}\\ \boldsymbol{{\phi}}_{y}^{st}\\ \boldsymbol{{\phi}}_{t}^{rot}\\ \boldsymbol{{\phi}}_{y}^{rot}\\ \boldsymbol{ i}^{st}\\ \boldsymbol{i }^{rot}\end{array}\right]-\left[\begin{array}{@{}c@{}}\boldsymbol{v}^{st}\\ \mathbf{0}\\ \boldsymbol{ 0}\\ \mathbf{0}\\ \boldsymbol{ 0}\\ \mathbf{0}\\ \boldsymbol{ 0}\\ \mathbf{0}\\ 0\\ 0\end{array}\right]=\mathbf{0}. & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$ (1) in a system of equations (SOE) given by Eq. (1), which is of the form $\begin{eqnarray}\displaystyle \boldsymbol{A}(\boldsymbol{x})\boldsymbol{x}-\boldsymbol{u}=\mathbf{0}. & & \displaystyle\end{eqnarray}$ (2) The convention is used that (i) matrices are denoted boldfaced and upper case, (ii) vectors are denoted boldfaced and lower case and (iii) scalars are denoted without boldface.

Fig. 1.

MEC model.

A ( x ) encompasses the system dynamics and is a function of the unknowns x . It is composed of various submatrices such as the winding distribution matrix N _abc. Further, resistance matrices are denoted using R and reluctance matrices using $\boldsymbol{{\mathcal{R}}}$ , where superscripts ^st and ^rot denote correspondence to the stator and rotor, respectively. Subscripts _t and _y denote correspondence to the teeth and yoke segments, respectively. Gauss’s law for fluxes is satisfied at the rotor and stator teeth edges using the air gap permeance submatrices $\boldsymbol{{\mathcal{P}}}_{st}^{rot},\boldsymbol{{\mathcal{P}}}_{rot}^{st}$ , and $\boldsymbol{{\mathcal{P}}}^{air}$ . I denotes an identity matrix, with its size given in subscript as either the number of stator teeth N_st or the number of rotor teeth N_rot. Lastly, the reference setting vector s e t _{1, N_st} ensures that the magnetic scalar potential of the first stator tooth is set to a reference value of 0. Further details of the MEC are provided in [13]. The nonlinear SOE is solved for x , which contains: stator tooth edge magnetic scalar potentials 𝜓^st, rotor tooth edge magnetic scalar potentials 𝜓^rot, stator tooth magnetic fluxes $\boldsymbol{{\phi}}_{t}^{st}$ , stator yoke magnetic fluxes $\boldsymbol{{\phi}}_{y}^{st}$ , rotor tooth magnetic fluxes $\boldsymbol{{\phi}}_{t}^{rot}$ , rotor yoke magnetic fluxes $\boldsymbol{{\phi}}_{y}^{rot}$ , stator currents i ^st and rotor bar currents i ^rot. Lastly, the input vector u contains the three phase-voltages v ^st, padded by zeros. The time-derivative in Eq. (1) is resolved using the trapezoidal rule [14].

2.2. Ferromagnetic saturation

The material-specific BH curve governs the magnetic permeability of the reluctances, where B is the magnetic flux density and H is the magnetic field strength. This material characteristic encompasses two nonlinearities, namely saturation and hysteresis. In this study, only saturation is considered. The used permeability stems from measured data of pure iron [15].

Figure 2 shows the permeability decrease due to saturation, where the material fully saturates around 1.7 T. The measured data was curve-fitted using a tanh function as $\begin{eqnarray}\displaystyle {\mu}(B)=-K_{{\mu}1}\tanh (K_{{\mu}2}|B|+K_{{\mu}3})+K_{{\mu}1}+{\mu}_{0}. & & \displaystyle\end{eqnarray}$ (3) Constants K_𝜇1 to K_𝜇3 are given in Table 1 and 𝜇₀ is the permeability of free space. The addition of 𝜇₀ ensures that the provided permeability remains larger than that of a vacuum. Equation (3) is used to express a derivative of the permeability with respect to the magnetic flux density. $\begin{eqnarray}\displaystyle {\displaystyle \frac{d{\mu}(B)}{dB}}=-\text{sgn}(B)K_{{\mu}1}K_{{\mu}2}\cosh (K_{{\mu}2}|B|+K_{{\mu}3})^{-2}. & & \displaystyle\end{eqnarray}$ (4)

The reluctances in A are expressed as a function of the magnetic flux. To illustrate this, submatrix $\boldsymbol{{\mathcal{R}}}_{t}^{st}$ is given by $\begin{eqnarray}\displaystyle \boldsymbol{{\mathcal{R}}}_{t}^{st}=(-{\delta}_{i,j}+{\delta}_{i,j-1})\text{diag}({\mathcal{R}}_{t,h}^{st})\quad \text{for }i,j,h\in [1,N_{st}], & & \displaystyle\end{eqnarray}$ (5) where 𝛿 is the Kronecker delta function [16] which uses i and j for indexing and ${\mathcal{R}}_{t,h}^{st}$ is the stator tooth reluctance [17] which is calculated as $\begin{eqnarray}\displaystyle {\mathcal{R}}_{t,h}^{st}={\displaystyle \frac{l_{\text{base}}}{S_{\text{base}}{\mu}\left(\frac{{\phi}_{t,h}^{st}}{S_{\text{base}}}\right)}}+{\displaystyle \frac{l_{\text{head}}}{S_{\text{head}}{\mu}_{\text{lin}}}}, & & \displaystyle\end{eqnarray}$ (6) with l_base and S_base corresponding to the length and the cross-section of the upper tooth part and l_head and S_head corresponding to the lower tooth part. 𝜇_lin is the linearised permeability, for which 𝜇(0) can be used. In this work, the tooth heads are assumed linear as they are unlikely to saturate due to their larger cross-sections. However, their nonlinearity can be included in an analogous manner to the first term of Eq. (6). Likewise, $\boldsymbol{{\mathcal{R}}}_{y}^{st}$ , $\boldsymbol{{\mathcal{R}}}_{t}^{rot}$ and $\boldsymbol{{\mathcal{R}}}_{y}^{rot}$ become dependent on the magnetic fluxes.

Table 1

Used solution parameters

Parameter	Usage	Value
K _𝜇1	Permeability function parameter	7.9 × 10⁻³
K _𝜇2	Permeability function parameter	4.4
K _𝜇3	Permeability function parameter	−6.1
K _𝛼1	Relaxation function parameter	0.90
K _𝛼2	Relaxation function parameter	−0.045
K _𝛼3	Relaxation function parameter	0.10
𝜖	Convergence threshold	1 × 10⁻²

Fig. 2.

Permeability of pure iron, curve fit and derivative of the permeability with respect to the magnetic flux density.

3. Nonlinear solution methods

As discussed in the previous section, A becomes dependent on x , yielding a nonlinear SOE which is expressed as $\begin{eqnarray}\displaystyle \boldsymbol{f}=\boldsymbol{A}\boldsymbol{x}-\boldsymbol{u}+\boldsymbol{c}=\mathbf{0}, & & \displaystyle\end{eqnarray}$ (7) where c comes from the trapezoidal rule applied to time discretisation.

3.1. Gauss–Seidel method (GSM)

GSM is an iterative method which is used to find the root of f . Let x _k denote the estimate of x at the kth iteration. Using a previous estimate of x , the estimate is updated as $\begin{eqnarray}\displaystyle \boldsymbol{x}_{k+1}^{\ast }=\boldsymbol{A}_{k}^{-1}(\boldsymbol{u}_{k}-\boldsymbol{c}_{k}), & & \displaystyle\end{eqnarray}$ (8) where e.g. A _k denotes the system matrix evaluated for x _k. Often, $\boldsymbol{x}_{k+1}^{\ast }$ as given by Eq. (8) is used directly to update the estimate. However, for substantial saturation levels, this does not converge. Instead, the produced estimate is only used in part, resulting in a relaxed version of GSM given by $\begin{eqnarray}\displaystyle \boldsymbol{x}_{k+1}=\boldsymbol{x}_{k}+{\alpha}(\boldsymbol{x}_{k+1}^{\ast }-\boldsymbol{x}_{k}), & & \displaystyle\end{eqnarray}$ (9) where 𝛼 is the relaxation factor. The relaxation factor is chosen as a function of the iteration number k such that relaxation increases as long as convergence remains unattained. $\begin{eqnarray}\displaystyle {\alpha}(k)=K_{{\alpha}1}\exp (K_{{\alpha}2}k)+K_{{\alpha}3}. & & \displaystyle\end{eqnarray}$ (10) The convergence is evaluated using the Frobenius norm of the residuals as $\begin{eqnarray}\displaystyle \Vert \boldsymbol{f}_{k}\Vert <{\epsilon}\Rightarrow \text{Converged}, & & \displaystyle\end{eqnarray}$ (11) where 𝜖 is a convergence threshold, given in Table 1. Since relaxation is applied, the change of x can be small even when the solution is not accurate. Therefore, the Frobenius norm of the residuals serves as a reliable measure of solution accuracy compared to the absolute or relative change of x .

3.2. Newton–Raphson method (NRM)

NRM uses knowledge of the derivatives of the SOE to perform a linearisation of the function in the current solution estimate. The root of this linearisation is used for $\boldsymbol{x}_{k+1}^{\ast }$ as $\begin{eqnarray}\displaystyle \boldsymbol{x}_{k+1}^{\ast }=\boldsymbol{x}_{k}-\boldsymbol{J}_{k}^{-1}\boldsymbol{f}_{k}, & & \displaystyle\end{eqnarray}$ (12) with $\boldsymbol{x}_{k+1}^{\ast }$ used in a damped manner as has been described in Eq. (9). The Jacobian J is a matrix which contains the partial derivatives of all equations to all unknowns as $\begin{eqnarray}\displaystyle \boldsymbol{J}=\left[\begin{array}{@{}cccc@{}}\frac{\partial f_{1}}{\partial x_{1}} & \frac{\partial f_{1}}{\partial x_{2}} & \cdots ∼ & \frac{\partial f_{1}}{\partial x_{m}}\\[3.60004pt] \frac{\partial f_{2}}{\partial x_{1}} & \frac{\partial f_{2}}{\partial x_{2}} & & \vdots \\[3.60004pt] \vdots & & \ddots & \vdots \\[3.60004pt] \frac{\partial f_{m}}{\partial x_{1}} & \cdots ∼ & \cdots ∼ & \frac{\partial f_{m}}{\partial x_{m}}\end{array}\right], & & \displaystyle\end{eqnarray}$ (13) where in general f_i is the ith equation of f and x_j is the j’th element of the m-length x . Using Eq. (7), the Jacobian is reduced to $\begin{eqnarray}\displaystyle \boldsymbol{J}=\left[\begin{array}{@{}cccc@{}}\frac{\partial \boldsymbol{A}_{1}\boldsymbol{x}}{\partial x_{1}} & \frac{\partial \boldsymbol{A}_{1}\boldsymbol{x}}{\partial x_{2}} & \cdots ∼ & \frac{\partial \boldsymbol{A}_{1}\boldsymbol{x}}{\partial x_{m}}\\[3.60004pt] \frac{\partial \boldsymbol{A}_{2}\boldsymbol{x}}{\partial x_{1}} & \frac{\partial \boldsymbol{A}_{2}\boldsymbol{x}}{\partial x_{2}} & & \vdots \\[3.60004pt] \vdots & & \ddots & \vdots \\[3.60004pt] \frac{\partial \boldsymbol{A}_{m}\boldsymbol{x}}{\partial x_{1}} & \cdots ∼ & \cdots ∼ & \frac{\partial \boldsymbol{A}_{m}\boldsymbol{x}}{\partial x_{m}}\end{array}\right]. & & \displaystyle\end{eqnarray}$ (14) It is easily shown that for all linear equations of f , (∂ A _i x ∕∂x_j) equals A _ij. Meanwhile, the nonlinear rows obtain an additional term which involves the partial derivative to x . This is exemplified by $\begin{eqnarray}\displaystyle \frac{\partial \boldsymbol{A}_{3}\boldsymbol{x}}{\partial \boldsymbol{{\phi}}_{t}^{st}}=\frac{\partial (\boldsymbol{{\mathcal{R}}}_{t}^{st}\boldsymbol{{\phi}}_{t}^{st}-\boldsymbol{{\mathcal{R}}}_{y}^{st}\boldsymbol{{\phi}}_{y}^{st}+\boldsymbol{N}_{abc}\boldsymbol{i}^{st})}{\partial \boldsymbol{{\phi}}_{t}^{st}}=\frac{\partial (\boldsymbol{{\mathcal{R}}}_{t}^{st}\boldsymbol{{\phi}}_{t}^{st})}{\partial \boldsymbol{{\phi}}_{t}^{st}}, & & \displaystyle\end{eqnarray}$ (15) where $\boldsymbol{{\phi}}_{t}^{st}$ and $\boldsymbol{{\phi}}_{y}^{st}$ are the stator teeth and yoke magnetic fluxes, $\boldsymbol{{\mathcal{R}}}_{t}^{st}$ and $\boldsymbol{{\mathcal{R}}}_{y}^{st}$ are the stator teeth and yoke reluctances, N _abc is the winding distribution matrix and i ^st are the stator currents. This is then evaluated for a general ${\phi}_{t,i}^{st}$ in $\boldsymbol{{\phi}}_{t}^{st}$ and ${\mathcal{R}}_{t,i}^{st}$ in $\boldsymbol{{\mathcal{R}}}_{t}^{st}$ as $\begin{eqnarray}\displaystyle \frac{\partial (-{\mathcal{R}}_{t,i}^{st}{\phi}_{t,i}^{st}+{\mathcal{R}}_{t,i+1}^{st}{\phi}_{t,i+1}^{st})}{\partial {\phi}_{t,i}^{st}}=-\frac{\partial ({\mathcal{R}}_{t,i}^{st}{\phi}_{t,i}^{st})}{\partial {\phi}_{t,i}^{st}}. & & \displaystyle\end{eqnarray}$ (16) Lastly, this is expanded using the chain rule as $\begin{eqnarray}\displaystyle -\frac{\partial {\mathcal{R}}_{t,i}^{st}}{\partial {\mu}_{t,i}^{st}}\frac{\partial {\mu}_{t,i}^{st}}{\partial {\phi}_{t,i}^{st}}{\phi}_{t,i}^{st}-{\mathcal{R}}_{t,i}^{st}. & & \displaystyle\end{eqnarray}$ (17) Where for a general reluctance ${\mathcal{R}}$ , its partial derivative with respect to 𝜇 equals $\begin{eqnarray}\displaystyle \frac{\partial {\mathcal{R}}}{\partial {\mu}}=\frac{\partial \frac{l}{S{\mu}}}{\partial {\mu}}=-\frac{l}{S{\mu}^{2}}, & & \displaystyle\end{eqnarray}$ (18) where l is the magnetic flux tube length and S is its cross-section. This process is repeated for all partial derivatives to the magnetic fluxes, which results in added terms involving the derivative of the permeability for the submatrices $\boldsymbol{{\mathcal{R}}}_{t}^{st}$ , $\boldsymbol{{\mathcal{R}}}_{y}^{st}$ , $\boldsymbol{{\mathcal{R}}}_{t}^{rot}$ and $\boldsymbol{{\mathcal{R}}}_{y}^{rot}$ .

3.3. Inverse Broyden’s method (IBM)

In Broyden’s method, the Jacobian is replaced by a difference quotient as $\begin{eqnarray}\displaystyle \boldsymbol{J}(\boldsymbol{x}_{k}-\boldsymbol{x}_{k-1})\approx \boldsymbol{f}_{k}-\boldsymbol{f}_{k-1}. & & \displaystyle\end{eqnarray}$ (19) For SOEs, Eq. (19) is underdetermined [18]. Broyden suggests updating J such that the Frobenius norm of its difference is minimised [18] as $\begin{eqnarray}\displaystyle \boldsymbol{J}_{k}=\boldsymbol{J}_{k-1}+{\displaystyle \frac{{\rm\Delta}\boldsymbol{f}_{k}-\boldsymbol{J}_{k-1}{\rm\Delta}\boldsymbol{x}_{k}}{\Vert {\rm\Delta}\boldsymbol{x}_{k}\Vert ^{2}}}{\rm\Delta}\boldsymbol{x}_{k}^{\intercal }, & & \displaystyle\end{eqnarray}$ (20) where Δ x _k = x _k − x _k−1 and Δ f _k = f _k − f _k−1. To improve the computational efficiency of Broyden’s method, the Sherman–Morrison formula [12] was applied to directly update the inverse of the Jacobian as $\begin{eqnarray}\displaystyle \boldsymbol{J}_{k}^{-1}=\boldsymbol{J}_{k-1}^{-1}+{\displaystyle \frac{{\rm\Delta}\boldsymbol{x}_{k}-\boldsymbol{J}_{k-1}^{-1}{\rm\Delta}\boldsymbol{f}_{k}}{{\rm\Delta}\boldsymbol{x}_{k}^{\intercal }\boldsymbol{J}_{k-1}^{-1}{\rm\Delta}\boldsymbol{f}_{k}}}{\rm\Delta}\boldsymbol{x}_{k}^{\intercal }\boldsymbol{J}_{k-1}^{-1}, & & \displaystyle\end{eqnarray}$ (21) which is called the Inverse Broyden’s Method (IBM). The Jacobian is calculated, inverted and stored using the Newton–Raphson method at the first iteration of each timestep. Thereafter, Eq. (21) is used to directly update the inverse, removing the need to further evaluate the exact Jacobian and to compute a matrix inverse. The estimate of x is then updated using Eqs (12) and (9).

4. Results and discussion

Magnetic saturation is examined using overloading and static rotor eccentricity. The simulation ran for 6000 time steps using a time step of 0.1 ms on a 2.6 GHz Intel Core i7-10750H CPU. The MEC configuration is based on a wye-connected, 230 V_RMS, 5 hp, 60 Hz motor. The used modelling parameters are given in Table 3. The model’s speed, active power, reactive power, power factor, output power, and power loss are evaluated during steady-state operation and compared with measurements in Table 4. The used nonlinear solution parameters are given in Table 1. The motor was loaded with 8 Nm of torque to induce a moderate global saturation and equipped with a rotor eccentricity of 50% of the air gap length in the positive x-direction, to induce a heavy local saturation. Depending on the figure either the last electrical period was used or the last 6 electrical periods were used to represent steady-state behaviour.

The saturation effects are validated using the magnetic flux density of the stator tooth in the direction of the eccentricity. The results of Fig. 3 agree with the material characteristic of Fig. 2, showing a reduction of the magnetic flux density to 1.74 T as the material’s permeability starts to approach that of free space. The figure also shows how the solutions obtained by the three solution methods are practically equal, as they are constrained by an equal tolerance on the norm of the residuals, see Eq. (11).

Fig. 3.

Magnetic flux density of the stator tooth base for h = 1 during last electrical period.

While the solution accuracy is equal, the difference between the methods lies in their performance. Figure 4 shows a considerable difference in the number of iterations over the last electrical period. GSM requires significantly more iterations when compared to NRM or IBM.

Fig. 4.

Number of iterations per time step for the three nonlinear solution methods during the last electrical period.

Figure 5 collects the required number of iterations for the entire steady state analysis. The observation is made that, on average, GSM requires the most iterations, followed by NRM, followed by IBM. Furthermore, analysing variance is valuable for control applications such as model predictive control where computational predictability is important. Here, the same order of performance is observed, i.e., the variance is largest for GSM, followed by NRM, and then IBM.

Fig. 5.

Number of iterations per time step for three nonlinear solution methods using a boxplot visualisation.

A performance comparison is given in Table 2. The results show the computation time per iteration to be the largest for NRM, followed by GSM, followed by IBM, although the differences per iteration are small. A general increase in the computation time per iteration is observed for nonlinear solution methods, which is attributed to the calculation of the norm of the residuals and other overhead. Furthermore, the computation time per iteration seems to be primordially determined by the calculation of inverse matrices and the Jacobian. IBM only explicitly calculates the inverse and the Jacobian once per time step, making it the most efficient solution method. GSM requires the calculation of the inverse of A , making it second. Lastly, NRM requires the calculation of the Jacobian and its inverse every iteration, making it the most costly. As IBM not only requires the smallest number of iterations but also the least computation time per iteration, it significantly outperforms NRM in terms of total computation time.

Fig. 6.

Convergence behaviour for a single time step.

Table 2

Performance comparison of the nonlinear solution methods

Method	Avg. #iterations	Computation time (s)	Time per iteration (ms)
No saturation	1	15.69	15.68
GSM	33.92	612.04	18.06
NRM	8.06	157.08	19.50
IBM	6.29	110.13	17.51

Lastly, Fig. 6 shows slow convergence in GSM caused by the relaxation factor, given by Eq. (9). For relaxation factors close to 0.5, the method becomes unstable and increases the error. Thereafter, as the relaxation factor drops below 0.45, the solution’s accuracy gradually increases again. This sensitivity to the chosen relaxation factor is an additional disadvantage. Further experimentation revealed NRM’s and IBM’s stability to be unaffected by relaxation. The error is generally reduced every iteration until the convergence threshold is reached. This happens consistently with fewer iterations for IBM than for NRM. A possible explanation for this is found in the derivative of the permeability as visualised in Fig. 2. For many magnetic flux density values, the derivative of the permeability is approximately zero. This results in an overshoot of NRM and is a known difficulty of the method in general. As IBM approximates NRM over two iterations, it is more likely to avoid a derivative of the permeability which is close to zero. Therefore, this method is less likely to overshoot the root and thus, on average, converges in fewer iterations.

Lastly, it is important to note that as saturation gets more severe, GSM is able to solve the system for the largest saturation levels, followed by NRM, and IBM. Therefore, for applications which require robust MEC simulations, multiple methods should be applied in a hybrid manner: e.g. IBM can be used by default, with a threshold on the number of iterations. If this threshold is exceeded, GSM can be invoked until it is able to solve the system in fewer iterations than a second threshold, whereafter the solution strategy reverts to IBM. A more extensive analysis of the convergence limitations is a topic for further research.

Table 3

Induction motor model parameters [3]

Model parameter	Value (unit)
Number of stator teeth	36
Stator inner diameter	104 mm
Stator outer diameter	170 mm
Stator tooth width	4.65 mm
Stator tooth head width	6.48 mm
Number of turns per slot	16.5
Stator tooth head thickness	0.2 mm
Stator phase resistance	406 mΩ
Pole pair number	2
Stator electrical connection	wye
Number of rotor teeth	28
Rotor inner diameter	35.0 mm
Rotor outer diameter	103 mm
Rotor slot depth	19.1 mm
Rotor tooth width	5.61 mm
Rotor tooth head width	10.6 mm
Rotor tooth head thickness	0.5 mm
Rotor end ring thickness	11.0 mm
Rotor axial length	120 mm
Rotor resistivity	3.37 × 10⁻⁸ Ωm
Rotor assembly inertia	0.01 kgm²
Rotor skew angle	5^°

Table 4

Model rotor speed and power comparison with reference to measurements [3]

	0.766 Nm load torque		24.7 Nm load torque
Quantity (unit)	Measurement	Simulation	Measurement	Simulation
Rotor speed (rpm)	1799	1799	1750	1753
Active power (W)	199.6	211.7	4966	4903.7
Reactive power (Var)	2609	2371	3426	3621
Power factor	0.076	0.089	0.823	0.804
Output power (W)	144.3	144.3	4563	4534
Power loss (W)	55.3	67.4	402.7	369.5

5. Conclusion

A MEC model of an induction motor is constructed, with saturation induced by overloading and rotor eccentricity. The nonlinearities, introduced by the magnetic flux-dependent reluctances in the stator and rotor, are elaborated. Three nonlinear solution methods are proposed for the comparison, namely: (i) the Gauss–Seidel method (GSM), (ii) the Newton–Raphson method (NRM) and (iii) inverse Broyden’s method (IBM). Of these methods, the first two have been explored extensively in literature. Solution accuracy is validated using the stator tooth magnetic flux density. Furthermore, a performance comparison revealed IBM to be the most performant solution method in terms of (i) the required number of iterations, (ii) the computation time per iteration and resultingly (iii) the total computation time. Compared to NRM, this is explained by the derivative of the permeability often being close to zero, which is a known weakness of NRM. Compared to GSM, this is explained by its sensitivity to the relaxation factor and its neglection of the system’s derivatives. However, for substantial saturation levels, the observation is made that GSM converges most, followed by NRM, and IBM. Therefore, the authors advise a hybrid implementation of multiple solution methods to obtain a flexible, efficient and robust solver. Future research should continue investigating the influences and limitations of the methods’ convergence.

Footnotes

Acknowledgements

This work was supported by the Research Foundation - Flanders (FWO) - application number 1S02523N and the TU Darmstadt research visit by CREATOR (DFG CRC TRR 361).

References

Naderi

Rostami

and Ramezannezhad

, Phase-to-phase fault detection method for synchronous reluctance machine using mec method, Electrical Engineering101(2) (2019), 575–586.

Ferreira

F.J.

Alberto

Silva

A.M.

and de Almeida

A.T.

, Saturation-related losses in induction motors for star and delta connection modes, in: 2020 International Conference on Electrical Machines (ICEM), Vol. 1, IEEE, 2020, pp. 1586–1593.

Sudhoff

S.D.

Kuhn

B.T.

Corzine

K.A.

and Branecky

B.T.

, Magnetic equivalent circuit modeling of induction motors, IEEE Transactions on Energy Conversion22(2) (2007), 259–270.

Tariq

A.R.

Nino-Baron

C.E.

and Strangas

E.G.

, Iron and magnet losses and torque calculation of interior permanent magnet synchronous machines using magnetic equivalent circuit, IEEE Transactions on Magnetics46(12) (2010), 4073–4080.

Tong

Wang

Dai

and Tang

, A quasi-three-dimensional magnetic equivalent circuit model of a double-sided axial flux permanent magnet machine considering local saturation, IEEE Transactions on Energy Conversion33(4) (2018), 2163–2173.

Bao

Gysen

and Lomonova

, Hybrid analytical modeling of saturated linear and rotary electrical machines: Integration of fourier modeling and magnetic equivalent circuits, IEEE Transactions on Magnetics54(11) (2018), 1–5.

Amrhein

and Krein

P.T.

, Induction machine modeling approach based on 3-d magnetic equivalent circuit framework, IEEE Transactions on Energy Conversion25(2) (2010), 339–347.

Tavana

N.R.

and Dinavahi

, Real-time nonlinear magnetic equivalent circuit model of induction machine on fpga for hardware-in-the-loop simulation, IEEE Transactions on Energy Conversion31(2) (2016), 520–530.

Ostovic

, Dynamics of Saturated Electric Machines, Springer Science & Business Media, 2012.

10.

Bash

M.L.

, Modeling and design of wound-rotor synchronous machines using mesh-based magnetic equivalent circuits, Ph.D. dissertation, Purdue University, 2010.

11.

Rabenstein

Dietz

Kremser

and Parspour

, Semi-analytical calculation of a laminated transverse flux machine, in: 2020 International Conference on Electrical Machines (ICEM), Vol. 1, IEEE, 2020, pp. 721–727.

12.

Maponi

, The solution of linear systems by using the sherman–morrison formula, Linear Algebra and its Applications420(2–3) (2007), 276–294.

13.

Desenfans

Gong

Vanoost

Gryllias

Boydens

and Pissoort

, The influence of the unbalanced magnetic pull on fault-induced rotor eccentricity in induction motors, Journal of Vibration and Control (2023).

14.

Rostami

Naderi

and Shiri

, Modeling and analysis of variable reluctance resolver using magnetic equivalent circuit, COMPEL-The International Journal for Computation and Mathematics in Electrical and Electronic Engineering40(4) (2021), 921–939.

15.

Boyer

H.E.

and Gall

T.L.

, Metals Handbook; Desk edition, ASM International, 1985.

16.

Kozen

and Timme

, Indefinite summation and the kronecker delta, 2007.

17.

Lannoo

Vanoost

Peuteman

Debruyne

De Gersem

and Pissoort

, Improved air gap permeance model to characterise the transient behaviour of electrical machines using magnetic equivalent circuit method, International Journal of Numerical Modelling: Electronic Networks, Devices and Fields33(5) (2020), e2749.

18.

Broyden

, On the discovery of the “good broyden” method, Mathematical Programming87 (2000), 209–213.