ON THE CALCULATION OF THE NONLINEAR MAGNETIC FIELD IN THE STATOR SLOT OF A SHOCK GENERATOR BY THE FINITE-DIFFERENCE METHOD

V. S. BAKLIN

(Presented by the scientific seminar of the Department of Electrical Machines and General Electrical Engineering)

The course of supertransient conditions in a shock generator, caused by a sudden short circuit, is substantially influenced by the leakage inductances of the rotor and stator circuits. Leakage inductances calculated without allowance for saturation differ significantly from the experimental ones obtained from tests of sudden short circuit at large current multiples [5]. Owing to the special design of shock generators, namely the presence of nonmagnetic rotor retaining rings and the shielding of the end portions of the stator winding by shields made of material with high electrical conductivity, it may be considered that the leakage inductances of the end portions of the rotor and stator circuits do not depend on saturation. Saturation has the principal effect on the slot-leakage inductance, since the transverse-slot leakage fluxes encounter on their path the steel regions of the stator and rotor teeth. Let us consider a method for calculating the nonlinear magnetic field in a stator slot.

The active length of the iron of the shock generator is assumed to be sufficient to neglect end effects and to consider that there is no variation of the field parameters along the \(z\)-axis (the longitudinal axis of the machine). This makes it possible to consider the leakage field in the slot in a two-dimensional space, assuming that the current-density vector \(\delta\) is directed along the machine axis.

Figure 1 shows the cross section of a symmetrical open slot. If one does not consider the leakage fluxes at the boundaries of the phase zones, assuming that the number of slots per pole and phase “\(q\)” is sufficiently large, then the field pattern will repeat periodically in the direction of the \(x\)-axis in each tooth pitch.

Since open symmetrical slots are considered, as the lines bounding the region under investigation we take the slot symmetry axis \(PE\), the tooth symmetry axis \(CD\), the middle of the air gap \(EC\), and the line \(PD\) bounding the bottom of the slot. To specify the boundary conditions it is necessary to establish the direction of the magnetic-field lines along the boundaries. The magnetic-field lines at the boundaries \(PE\) and \(CD\) will have only one normal component, since \(PE\) and \(CD\) are lines of symmetry. The magnetic permeability of the yoke, in comparison with the magnetic permeability of the teeth, is taken to be infinite. Consequently, the magnetic-field lines must be perpendicular to the boundary line \(PD\).

In the case of an ideal short circuit by the mutual-induction flux between—

Fig. 1. For the calculation of the magnetic field in the slot

the rotor and stator contours can be neglected, and the leakage fields of the stator and rotor can be regarded as unconnected with one another. When the m.m.f.’s of the stator and rotor contours are equal, the leakage fluxes of each of them reach approximately the middle of the air gap (up to line \(EC\), and do not cross it). We take line \(EC\) to be the line of zero potential. The bars in the stator slot of a shock generator, with rare exceptions, consist of elementary conductors transposed among themselves; consequently, the influence of eddy currents may be neglected. We neglect the phenomenon of hysteresis in the tooth, assuming that the magnetization curve of the steel is determined uniquely.

To calculate the nonlinear magnetic field in the slot, we shall use the finite-difference method, which belongs to the class of numerical methods. The essence of this method consists in approximating the field equation by a system of finite-difference equations, and the calculation of the magnetic field is reduced to solving this system. General questions of the finite-difference method are set forth in [1, 2], and as applied to the calculation of nonlinear magnetic fields in electrical machines in [4, 7, 8, 9].

Equations in finite differences can be obtained by using the law of total current in integral or differential form [4]. If in the region under investigation there are boundaries of separation of media where the magnetic permeability \(\mu\) and the current density \(\delta\) change discontinuously, then it is necessary to use the law of total current in integral form, which is conveniently reduced to a form containing, as the unknown function, the vector potential of the magnetic field \(A\).

\[ \oint \frac{1}{\mu}\operatorname{rot}\vec{A}\cdot d\vec{l} = \int_s \vec{\delta}\cdot d\vec{s}. \tag{1} \]

The calculation equation for the vector potential, obtained on the basis of

expressions (1), as a result of traversing the contour \(1, 2, 3, 4, 5, 6, 7, 8, 1\), Fig. 2, for each internal grid node \(i, k\), has the form

\[ A_{i,k}=\frac{1}{M_{i,k}}\left(N_{i,k}\cdot A_{i,k+1}+S_{i,k}\cdot A_{i,k-1}+O_{i,k}\cdot A_{i+1,k}+\right. \]
\[ \left.+\,W_{i,k}\cdot A_{i-1,k}+D_{i,k}\right); \tag{2} \]

\[ M_{i,k}=N_{i,k}+S_{i,k}+O_{i,k}+W_{i,k}; \tag{3} \]

\[ N_{i,k}=\frac{1}{2\cdot q_k}\cdot \left(\frac{p_i}{\mu_{i,k}}+\frac{p_{i-1}}{\mu_{i-1,k}}\right), \quad S_{i,k}=N_{i,k-1}; \tag{4} \]

\[ O_{i,k}=\frac{1}{2\cdot p_i}\cdot \left(\frac{q_k}{\mu_{i,k}}+\frac{q_{k-1}}{\mu_{i,k-1}}\right), \quad W_{i,k}=O_{i-1,k}; \tag{5} \]

\[ D_{i,k}=\int \delta\cdot ds= \frac{1}{4}\left(\delta_{i,k}\cdot p_i\cdot q_k+\delta_{i-1,k}\cdot p_{i-1}\cdot q_k+\delta_{i-1,k-1}\times\right. \]
\[ \left.\times\,p_{i-1}\cdot q_{k-1}+\delta_{i,k-1}\cdot p_i\cdot q_{k-1}\right), \tag{6} \]

where \(p_i, q_k\) are the dimensions of the grid cells,

\(\mu_{i,k}, \delta_{i,k}\) are the magnetic permeability and current density in the grid cells,

\(A_{i,k}\) are the values of the vector potential at the grid nodes.

The value of the current density \(\delta\) and the magnetic permeability \(\mu\) within each grid cell is assumed to be constant. For this it is necessary that the grid lines coincide with the interface between media, where the magnetic permeability and current density change discontinuously.

Fig. 2. Arrangement of the calculated quantities and their designations on a rectangular grid

In calculating nonlinear magnetic fields it is necessary to know the modulus of the magnetic induction vector, which can be found in the following way:

\[ B_{x\,i,k}=\frac{1}{2\cdot q_k}\cdot \left(A_{i,k+1}+A_{i+1,k+1}-A_{i,k}-A_{i+1,k}\right); \tag{7} \]

\[ B_{y\,i,k}=-\frac{1}{2\cdot p_i}\cdot (A_{i,k+1}+A_{i,k}-A_{i+1,k+1}-A_{i+1,k}); \tag{8} \]

\[ B_{i,k}=\sqrt{B_{x\,i,k}^{\,2}+B_{y\,i,k}^{\,2}}. \tag{9} \]

The magnitude of the magnetic induction vector \(B_{i,k}\) at the center of the cell is used to determine the magnetic permeability \(\mu_{i,k}\) from the iron magnetization curve.

When calculating magnetic fields with allowance for saturation, it becomes necessary to approximate the iron magnetization curve. In view of the complicated dependence \(B=f(H)\), it is practically impossible to approximate the entire iron magnetization curve by a single simple analytic expression. Most often a piecewise approximation is used, in which each segment of the magnetization curve can be approximated by simple algebraic or transcendental functions. If the number of segments is large, then one may use the method of piecewise-linear approximation [10], which is quite easily implemented by means of computers. The entire magnetization curve up to complete saturation (up to the induction \(B_s\)) is divided into segments with step \(\Delta B\). The magnetic permeability of iron for \(B<B_s\) is determined from the expression

\[ \mu_{i,k}=\frac{B_{i,k}}{(B_{i,k}-B_j)\cdot m_j+H_j}, \tag{10} \]

where

\[ m_j=\frac{H_{j+1}-H_j}{\Delta B} \]

for each segment \(j\),

and the segment number is

\[ j=\frac{B_{i,k}}{\Delta B}+1. \]

For \(B\ge B_s\)

\[ \mu_{i,k}=\mu_0\cdot \frac{B_{i,k}}{B_{i,k}-B_{\mathrm{kp}}}. \tag{11} \]

A finite-difference approximation of equation (1) gives a system of nonlinear algebraic equations, the number of which is equal to the number of internal nodes of the grid. The coefficients of the unknowns in the system of nonlinear algebraic equations form a sparse band matrix. For solving large systems of algebraic equations with a sparse matrix, iterative methods are very convenient. Among iterative methods, the Gauss–Seidel method is used most widely. In connection with elliptic difference equations it is often called the Liebmann method [1, 2, 8].

If special acceleration measures are not applied, the iteration method is characterized by very slow convergence. In calculating magnetic fields by the finite-difference method, two convergence-acceleration methods are most widely used, which involve:

1) using the total-current law written in integral form [8, 9];

2) using incomplete upper relaxation [1, 2, 3].

In the first case, the total-current law written in integral form,

\[ \oint \bar{H}\cdot d\bar{l}=\int \bar{\delta}\cdot d\bar{s} \tag{12} \]

is used as a means of checking and accelerating convergence. After the first iterations, equation (12) is usually not satisfied exactly.

\[ C=\frac{\displaystyle\int_s \bar{\delta}\cdot d\bar{s}}{\displaystyle\oint \bar{H}\cdot d\bar{l}}\ne 1. \tag{13} \]

If the integration contour \(G\) is laid out in such a way that it coincides with some field line, where \(A_G=\mathrm{const}\), then the coefficient “\(C\)” can be used to accelerate the convergence of the iterative process. In calculating the contour integral \(\oint \bar{H}\cdot d\bar{l}\), the value of the minimum potential \(A_{\min}\) along the integration contour is determined. Then all potentials with \(|A|>|A_{\min}|\) are corrected according to the following rule:

\[ A'=(A-A_{\min})\cdot C+A_{\min}. \tag{14} \]

If possible, the integration path is chosen so that \(A_{\min}=0\), which is often determined by the boundary conditions. In this case expression (14) takes the form

\[ A'=A\cdot C, \tag{15} \]

where \(A'\) is the corrected value of the vector potential.

In addition to the coefficient “\(C\),” the value of the vector potential \(\Delta A\), obtained from the expression

\[ \Delta A= \frac{\displaystyle\int_s \bar{\delta}\cdot d\bar{s}-\displaystyle\oint \bar{H}\cdot d\bar{l}}{L}, \tag{16} \]

is also used to accelerate convergence.

In the case of rectangular grid cells, the values \(L\), \(\oint \bar{H}\cdot d\bar{l}\), and \(\int_s \bar{\delta}\cdot d\bar{s}\) can be determined as follows:

\[ L=\sum_{G_x}\frac{p_i}{q_k\cdot \mu_{i,k}} +\sum_{G_y}\frac{q_k}{p_i\cdot \mu_{i,k}}; \tag{17} \]

\[ \int_s \bar{\delta}\cdot d\bar{s}=I=\sum \delta_{i,k}\cdot p_i\cdot q_k; \tag{18} \]

\[ \oint \bar{H}\cdot d\bar{l} = \sum_{G_x} \left( \frac{A_{i,k+1}+A_{i+1,k+1}-A_{i,k}-A_{i+1,k}} {2\cdot q_k\cdot \mu_{i,k}} \right)\cdot(\pm p_i) + \]
\[ + \sum_{G_y} \left( \frac{A_{i,k}+A_{i,k+1}-A_{i+1,k}-A_{i+1,k+1}} {2\cdot p_i\cdot \mu_{i,k}} \right)\cdot(\pm q_k). \tag{19} \]

In expression (19), the sign \((+)\) is used in the case when the integration path coincides with the direction of the unit vectors \(i\), \(j\) (Fig. 1). The correction of the vector potential within the integration path for expression (16) is performed according to the formula

\[ A'=A+\Delta A. \tag{20} \]

The essence of the method for accelerating convergence, which uses incomplete over-relaxation, consists in the way the error “\(r\)” is suppressed for the unknowns

\[ r=A_{i,k}^{n}-A_{i,k}^{(n-1)}, \tag{21} \]

where \(n\) is the iteration number.

Depending on the method of suppressing the error “\(r\),” relaxation may be complete or incomplete. The type of relaxation is determined by the value of the coefficient \(\omega\) in the expression

\[ A_{i,k}^{(n+1)}=A_{i,k}^{(n-1)}\cdot(1-\omega)+\omega A_{i,k}^{n}, \tag{22} \]

where \(A_{i,k}^{n}\) is determined from the equation

\[ A_{i,k}^{n}=\frac{1}{M_{i,k}}\cdot\left(N_{i,k}\cdot A_{i,k+1}^{n}+S_{i,k}\cdot A_{i,k-1}^{n}+W_{i,k}\cdot A_{i-1,k}^{(n+1)}+\right. \]
\[ \left.+\,O_{i,k}\cdot A_{i+1,k}^{n}+D_{i,k}\right), \tag{23} \]

if the Gauss--Seidel method is used. According to [3], if \(\omega=1\), then complete minimization of the error function “\(r\),” or complete upper relaxation, takes place. Incomplete relaxation is called lower relaxation if \(0<\omega<1\), and upper relaxation if \(1<\omega<2\). To solve a system of finite-difference equations by the Gauss--Seidel method it is advantageous to apply incomplete upper relaxation, which is associated with the tridiagonal representation of the matrix of coefficients of the vector potentials, possessing property A [2]. The rate of convergence of the iteration process depends strongly on how close the value of the coefficient of incomplete upper relaxation is to the optimal one. If the region under study has a complex boundary configuration, then there are no analytical methods for determining the optimal value of the coefficient of incomplete upper relaxation, and its value is established experimentally.

Sometimes, in the course of the calculation, automatic optimization of the coefficient \(\omega\) is applied [4, 7]. For this, in addition to carrying out auxiliary calculations, it is necessary to store the values of the vector potentials for three successive iterations, which substantially overloads the computer’s random-access memory.

In order for the iteration process to proceed stably, in nonlinear problems the magnetic permeability is usually corrected:

\[ \mu_{i,k}^{(n+1)}=\mu_{i,k}^{n}\cdot(1-\alpha)+\alpha\cdot\mu_{i,k}, \tag{24} \]

where \(\mu_{i,k}\) is determined from expression (11) or (12).

In each specific case the values of the coefficient \(\alpha\) are established experimentally. According to [7, 8], the iteration process proceeds stably if \(\alpha\leq 0.1\). In [9] the coefficient \(\alpha\) varies within the range \(0.25\div 0.0625\).

To calculate the magnetic field, the region under study is covered by a grid which, in order to save the computer’s random-access memory, is often made nonuniform.

Figure 3 shows the region of one tooth pitch covered by a grid. The grid has 378 computational nodes. The dimensions of the slot and tooth are as follows:

1) slot width \(b_{\mathrm{p}}=5\ \mathrm{cm}\);

2) tooth width \(b_{z}=5\ \mathrm{cm}\);

3) slot height \(h_{\mathrm{p}}=10\ \mathrm{cm}\);

4) height of the wedge together with the insulation \(h_{k}=2\ \mathrm{cm}\);

5) magnitude of the air gap up to line ES: \(1.5\ \mathrm{cm}\).

The stator is made of electrical steel E41, for which \(B_{\mathrm{kr}}=1.914\ \mathrm{T}\), \(\Delta B=0.1\ \mathrm{T}\), \(B_{s}=2.1\ \mathrm{T}\) were adopted. At inductions \(B\leq 0.1\ \mathrm{T}\), the magnetic permeability of the steel was taken equal to \(0.157\cdot 10^{-3}\ \mathrm{H/m}\).

By expression (23) only the values of the vector potentials of the internal cells of the grid are calculated, but not those on the boundary of the region under study. To calculate the vector potentials at the grid nodes lying on the boundary of the region under study, it is necessary to impose boundary conditions on the lines RE, SD, ES, RD. The fulfillment of the boundary conditions on the lines RE, SD, RD is ensured by introducing artificial boundaries, which, according to Fig. 3, are the grid lines \(0,k\) for RE, \(15,k\) for SD, and \(i,0\) for RD. After each iteration the vector potentials on the artificial boundaries are assigned the values

\[ A_{0,k}=A_{2,k};\qquad A_{15,k}=A_{13,k};\qquad A_{i,0}=A_{i,2}. \]

Fig. 3. The investigated region divided into cells of a rectangular grid.

The line \(EC\) is taken as the line of zero potential; consequently,

\[ A_{ec}=A_{1,29}=0. \]

The magnetic-field calculation was carried out on an “M-222” computer. The program was written in ALFA-language [6].

To accelerate convergence of the iterative process, both the total-current law written in integral form and incomplete upper relaxation were used. The system of nonlinear finite-difference equations was solved by the Gauss-Seidel method. The calculation of the vector potentials began from node \((1,1)\) in Fig. 3 and was carried out by columns. The solution was considered to have been reached if, for all internal nodes of the grid, the condition

\[ \Delta=\operatorname{mod}\left(A_{i,k}^{(n+1)}-A_{i,k}^{n}\right)\leq 10^{-3},\ 10^{-4},\ 10^{-5}. \tag{25} \]

was satisfied.

If, to accelerate convergence, the total-current law in integral form was used, then the coefficient \(C\) was calculated from expression (13). Integration was carried out along the contour \(G_x\) (Fig. 3). When condition (25) was fulfilled, the program provided for the output of the following results:

1) the values of the vector potential at all grid nodes;
2) the values of the magnetic permeability in the cells of the tooth region;
3) the number of iterations “IT”;
4) the value of the coefficient “\(C\)”;
5) the value of the total slot conductance coefficient \(\lambda_{\mathrm{tot}}\).

As initial values, zero values of the vector potential were taken and, for the cells in the tooth region, values of the magnetic permeability \(\mu=200\mu_0\). The studies showed that applying the total-current law in integral form to accelerate convergence of the iterative process over the entire region bounded by the integration contour leads to instability of the iterative process. If, however, the total-current law is used only to correct the vector potentials within the integration contour, then the iteration process proceeds relatively stably, although the computation time increases substantially compared with the method of accelerating convergence of the iterative process that uses incomplete upper relaxation. Some results of the calculation of the magnetic field for a current density \(\delta=150\ \mathrm{A/mm^2}\) are presented in Table 1. The values in the numerator refer

Table 1

to the method of accelerating convergence of the iterative process that uses the total-current law in integral form, and those in the denominator to the method using incomplete upper relaxation. According to Table 1, with an accuracy acceptable for engineering calculations, the solution may be considered reached if

\[ \Delta=\operatorname{mod}\left(A_{i,k}^{(n+1)}-A_{i,k}^{n}\right)\leq 10^{-5} \]

or

\[ \operatorname{mod}(1-C)\leq 10^{-3}. \]

At the same time, about 50% of the machine time is spent on reducing the modulus of the error \(r\) from \(10^{-4}\) to \(10^{-5}\). From the standpoint of machine-time expenditure, the most advantageous in this case is the use of incomplete upper relaxation. The slow convergence of the iterative process, when the law of full current is used in integral form, can be explained by the fact that the coefficient \(C\) rapidly tends to unity, and in the final stages, where the modulus of \(r\) decreases from \(10^{-3}\) to \(10^{-5}\) and the main machine time is spent, it has no substantial effect on the rate of convergence of the iterative process. Solving the given problem with the use of incomplete upper relaxation required 22 minutes of machine time, of which about 10 minutes were spent on program translation. The desire to reduce the computation time led to the choice, for the same investigated region, of a mesh with 198 internal nodes. A comparison of the calculation results for the total slot conductance coefficient \(\lambda_{\mathrm{tot}}\) for the two meshes is presented in Table 2, where the data under number 1 refer to the mesh with 378 nodes, and the data under number 2 correspond to the mesh with 198 nodes. For both meshes \(\omega=1.75\), \(a=0.08\). The difference in the calculation results for \(\lambda_{\mathrm{tot}}\) does not exceed 5%, whereas the machine-time expenditures were, for the mesh with 378 nodes for all current densities \(\delta\) presented in Table 2, about

Table 2

90 minutes, and for the mesh with 198 nodes, 22 minutes. For the mesh with 198 internal nodes, the influence of the value of the coefficient of incomplete upper relaxation on the convergence rate of the iterative process was investigated. It was observed that the optimal value of the coefficient of incomplete upper relaxation, at which the maximum convergence rate of the iterative process is achieved, depends on the ratio of the regions occupied by iron and free of iron, on the degree of saturation of the iron, on the adjustment of the coefficient \(a\), on the geometrical dimensions of the mesh cells, etc. Thus, for example, by reducing the region occupied by iron, one can somewhat increase the value of the coefficient of incomplete upper relaxation. This is clearly seen from Tables 3 and 4. The use of the coefficient of incomplete upper relaxation \(\omega=1.85\) instead of \(\omega=1.75\), and adjustment of the coefficient \(a\) during the iterations, reduce the expenditure of machine time by 30%. In this case the tooth width is \(b_z=2.5\) cm. Similar investigations for \(b_z=10\) cm showed that the iteration process diverges at \(\omega=1.85\). For the mesh with 198 nodes, the method was also investigated

Table 3

Table 4

variable directions [2]. No evident advantages of the latter method over the method using incomplete upper relaxation were observed, which is apparently connected with the relatively small number of grid nodes, as well as with the nonlinearity of the problem.

Conclusions

If there is no need for a graphical construction of the magnetic-field pattern, then it is quite sufficient, in the region under study, occupying half of the stator tooth pitch, to restrict oneself to 198 internal nodes. In this case, the Gauss–Seidel iteration method using incomplete upper relaxation is implemented most simply on a computer.

References

1. K. Binns, P. Lawrenson. Analysis and Calculation of Electric and Magnetic Fields. Moscow, “Energiya,” 1970.

2. V. Vazov, Dzh. Forsayt. Difference Methods for Solving Partial Differential Equations. Moscow, IL, 1963.

3. D. K. Fadeev, V. N. Fadeeva. Computational Methods of Linear Algebra. Moscow, Fizmatgiz, 1963.

4. Ya. A. Novik. Numerical methods for calculating the magnetic field of electrical machines taking saturation into account. In: Contactless Electrical Machines, 11, Riga, “Zinatne,” 1972.

5. L. F. Bukreev, N. V. Shilin. The influence of saturation on the inductive resistances of shock generators of the GI-100-2 type, Elektrotekhnika, 1970, 11.

6. A. P. Ershov, G. I. Kozhukhin, I. V. Pottosin. Guide to the Use of the ALFA System. Novosibirsk, “Nauka,” Siberian Branch, 1968.

7. K. Reichert. Ein numerisches Verfahren zur Berechnung magnetische Felder, insbesondere in Anordnungen mit Permanentmagneten. Arch. Elektrotech., 52, 1968, 3.

8. K. Reichert. Die Numerische Berechnung der Nutstreuung und der Wirbelströme Nutenleitern elektrischen Maschinens. Arch. Elektrotech., 53, 1969, 2.

9. H. Braess, H. Weh, E. A. Erdelyi. Numerische Berechnung magnetischer Felder und Kräfte. Arch. Elektrotech., 52, 1969, 5.

10. Frederick C., Erdelyi E. A., Hopkins R. E. Representation of the Magnetization Characteristic of DC Machines for Computer Use. — IEE Trans. Power Apparatus and Systems, Pas—87, 1968, 3.