Evaluation of High-order Harmonic Currents of Induction
Motors Using A Variant OF Crank-Nicholson’s Method
S. L. Ho W. N. Fu
Department of Electrical Engineering
The Hong Kong Polytechnic University, Hong Kong
Email:
Abstract
A numerical model for evaluating high-order harmonic currents of skewed rotor induction motors is described. The model is based on the time-stepping finite element method. When the applied voltage is known, the true waveform of the stator current can be computed directly. In discretizing the time variable in the proposed algorithm, a variant of Crank-Nicholson’s method is presented, and the proposed method has the merit of being as simple as the Backward Euler’s method with a better accuracy. The computed results are verified by the test results.
1. INTRODUCTION
Induction motors are one of the main loads in power networks. Because of the existence of slots in stator and rotor, the non-sinusoidal arrangement of magnetomotive force in space and the saturation of iron cores, it is commonly observed that the stator current may contain many high-order harmonic components. These harmonics will certainly affect the quality of power networks. On the other hand, the study of high-order harmonics is also the first step in studding stray losses, additional torques, vibration and noise. Attempt to evaluate correctly the magnitudes of high-order harmonic currents is therefore an important challenge for electrical machine designers.
The traditional method based on the magnetic circuit concept cannot solve such problem precisely. The normal commercial finite element software can produce the distribution of electromagnetic field only if the stator current is known. The application of time stepping finite element method which is coupled with external circuit equations and the electromagnetic equations, and which moves the finite element mesh in sympathy with the rotation of the rotor, is indeed an excellent method in estimating the non-sinusoidal waveforms of the stator current [1-3]. However, this method requires a large amount of computing time. In this paper a variant of the Crank-Nicholson’s (C-N) method is presented. This approach has the merits of having a simple algorithm compared with the standard C-N method and it requires very little computing time compared with that of the Backward Euler’s method for the same allowable error.
The method is applied to study the on load operation of an 11 kW, cage, skewed induction motor. The stator current waveforms at different skewed angles of rotor bars are computed and compared.
2. MATHEMATICAL MODEL OF
INDUCTION MOTORS
The motor in the axial direction is considered as composing of multi-slices [3-4] (Fig. 1). In each slice the magnetic vector potential has an axial component only. The magnetic fields are present in planes normal to the machine axis. Hence the characteristic of the electromagnetic field of each slice is two-dimensional. The electrical relationships between the various slices are based on the principle that the current flowing in the bars of one slice is the same as that which flows in the bars of every other slices.
(a) A skewed rotor (b) A multi-slice model
Fig. 1 Model for the representation of skew
(only one quarter of rotor is shown)
The stator outside circuit and the effect of the rotor end rings are considered by coupling the electrical circuits into the FEM equations. The leakage inductances at the end of the stator winding and at the end ring of the rotor cage are obtained by analytical methods.
The Maxwell’s equations applied to the domain under investigation will give rise to the following equation:
(1)
where A is the axial component of the magnetic vector potential; n is the reluctivity of the material and j is the total current density and includes the eddy current.
Because the iron cores are laminated, the eddy currents in the iron cores can be neglected in the mathematical model. Therefore, in the iron domains and in the air-gap, .
In the rotor conductor domain, the total current density j(k) in the conductor of the kth slice is:
(2)
where l(k) is the axial length of the kth slice, u(k) is the potential difference between the two ends of a conductor in the kth slice and s is the conductivity of the material.
Integrating the current density j(k) over the cross-section of the conductor gives:
(3)
where i(k) is the total current in the conductor. S is the cross-sectional area of the rotor bar.
By putting l(k) in equation (3) to the left hand side, and letting k=1,2,...,M (M is the number of slices) and to group all these M equations together, one has
(4)
Replace i(m) in equation (4) with i, one obtains:
(5)
where is the axial length of the bar and is the potential differences between the two ends of the bar.
Because the items on the right hand side of equation (3) and (5) are equal, one has
(6)
Substituting (6) into (2), and noting j(k) , A(k) and W(k) can be written simply as j, A and W, one obtains:
(7)
Another relationship between i and u can be obtained by the circuit equations for the rotor end windings [2]. The unknown potential difference u can then be eliminated at first. Equations (1), (7), (5) and the end winding circuit equations will give rise to the basic formulas in the rotor bar domain.
In the stator conductor domain, because the stator winding consists of fine wires, the current density j in the stator conductor is assumed uniform.
Assuming the stator winding of one phase consists of w turns in series, the cross-sectional area of one phase conductors at one side of the coils is W+, the area at the other side of the coils is W- and the total cross-sectional area of one turn at one side is equal to S. The circuit equations are:
(8)
(9)
where Vs is the applied voltage; i is the phase current; Rs and Ls are resistance and inductance at the end of winding, respectively; u is the total potential difference of each element in the stator winding; R is the resistance of each element in the stator. The total induced electromotive force in one phase is therefore [5]:
(10)
Substituting equation (10) into (9), and noting that , the total conductor area of one phase is and , one obtains
(11)
Substituting equation (8) into (11), one has
(12)
Equation (1) and (12) constitute the basic formulas governing the stator conductor domain.
By using finite element formulation and coupling the electromagnetic field equations together with the electrical circuit equations, one obtains the following large nonlinear system of equations:
(13)
where the unknowns [A] and [i] are the magnetic vector potentials and currents, respectively, that are required to be evaluated. [K], [C], [Q], [R] are the coefficient matrices and [P] is the vector associated with input voltages.
3. BACKWARD EULER’S METHOD AND THE
CRANK-NICHOLSON’S (C-N) METHOD
Equation (13) can be simply written as:
(14)
where X is an unknown vector that is required to be solved and that include the magnetic vector potentials [A] and the currents vector [i]; D and B are coefficient matrices with sparse characteristics.
Because the Backward Euler’s method (with first order precision) and the Crank-Nicholson’s (C-N) method (with second order precision) are A-stable [6][7], they are being used frequently to discretize the time variable in finite element problems [1][8-10].
In the Backward Euler’s method,
(15)
Notice that is determined by (14), Xk can be achieved by the solution of the following algebraic equations:
(16)
where the step size .
In the C-N method,
(17)
where is determined by
(18)
and is determined by
(19)
Because the global discretization error of the Backward Euler’s method satisfies , and that of the C-N method is , it is apparent that the C-N method can use larger step size in Dt. In achieving the same accuracy, it is worth pointing out that although the computation time of the C-N method is a little more than that required by the Backward Euler’s method at each step, the total computing time is still less for the latter. However, the difficulty is, in the combined equation incorporating (17), (18) and (19), Dk-1, Bk-1 and Dk, Bk exist all at the same time. Note that because the mesh of the rotor part is being rotated, the relationship of nodes and elements in the mesh at the (k-1)th step is different with that at the kth step. Since Dk-1 and Bk-1 are associated with the mesh at the (k-1)th step while Dk and Bk belong to the mesh at the kth step, the sparse conditions of Dk-1 and Bk-1 are different from that of Dk and Bk. The coefficient matrixes will also be stored in compressed style, and this is very essential in saving computer RAM and CPU time. Therefore, the algorithm required to solve the combined equation (17), (18) and (19) will become very complicated.
4. VARIANT C-N METHOD
To simplify the problem, the authors will describe a variant form of the C-N method. Firstly, the mesh will be generated according to the following method so that the matrix B can be kept constant at the different time steps.
Firstly, the domain to be investigated is divided into two parts at the middle of the air gap to produce the stator and the rotor domains. In other words, each of domains includes one half of the corresponding air gap in its vicinity. The stator meshes and the rotor meshes are then generated separately. The nodes at the inner-most of the stator and at the outer-most of the rotor are arranged in equal distance. The unknowns at the inner-most nodes of the stator and at the outer-most nodes of the rotor are connected according to “periodic boundary conditions”. That is to say, when the rotor is rotated, the shape of the rotor meshes will not be changed, it is only the co-ordinates and the periodic boundary conditions which are to be changed. Thus the stator meshes and the rotor meshes need to be generated once only. The node numbers at all nodes are kept unchanged. The other important merit of this mesh generation method is, because there is no eddy current in the air gap and B is the coefficient matrix for , the change of the boundary conditions will not change the matrix B when the rotor meshes are rotated. That means:
(20)
Multiplying B at both sides of (17), gives:
(21)
Substituting and from (18) and (19) into equation (21) gives:
(22)
where the vector can be obtained directly from the k-1 step.
Comparing (16) and (22), it is clear that the algorithm of the C-N method becomes as simple as the Backward Euler’s method. It is very easy to modify the Backward Euler’s method to implement the C-N method in the software.
5. NUMERICAL RESULTS
The presented method has been used to simulate the steady-state on load operation of an 11 kW, skewed rotor cage induction motor. The details of the motor are shown at table 1.
Table 1 Main parameters of the test motor
Rated power (kW) 11Rated voltage (v) 380
Rated current (A) 22.61
Rated frequency (Hz) 50
Connection D
Number of phases 3
Number of pole pairs 2
Number of stator slots 48
Number of rotor slots 44
Stator diameter (mm) 240
Rotor diameter (mm) 157
Air-gap length (mm) 0.5
Core length (mm) 165
Skew 1.3 stator slot pitch
The program is run on a Pentium / 90 MHz personal computer. The size of each time step is 0.078 ms. The slice number is 5. The domain being investigated is one pair of poles (Fig. 2). Each slice contains 3232 nodes and 5304 elements. The total number of unknowns at each step is 16185. The computation of each step needs a CPU time about 190 s.
Fig. 2 A FEM mesh at the one slice of 11 kW
skewed motor after rotation of 198.41°
Noted that a high-order harmonic index H is introduced to show the magnitude of high-order harmonic components:
% (23)
where I(k) is the magnitude of the kth harmonic current waveform expanded in Fourier series.
One of the computed flux distribution is shown in Fig. 3. The computed stator phase current is shown in Fig. 4, with H=7.0 %. The stator phase current waveform obtained by experiment is shown in Fig. 5, its H is 6.8 %. It is obvious that the results obtained by using the algorithm being developed has a very good correlation with the test data.
Fig. 3 The flux distribution at the one slice of 11 kW
skewed motor after rotation of 198.41°
Fig. 4 Computed stator phase current waveform
(Skew of 1.3 stator slot pitch)
Fig. 5 Measured stator phase current waveform
(Skew of 1.3 stator slot pitch)
The program is also used to compute the stator phase currents when the rotor bars are skewed at different angles. The results are shown in Fig. 6, Fig. 7, Fig. 8. and table 2. It can be seen that when the skewed angle is about one stator slot pitch, the high-order harmonics can be reduced.
Fig. 6 Computed stator phase current waveform
(Non-skewed Motor)
Fig. 7 Computed stator phase current waveform
(Skew of 0.7 stator slot pitch)
Fig. 8 Computed stator phase current waveform
(Skew of 1 stator slot pitch)
Table 2 Computed high-order harmonic index H
when the rotor bars are in different
skewed angles
Skew 0 stator slot pitch H=10.3 %Skew 0.7 stator slot pitch H= 7.4 %
Skew 1 stator slot pitch H= 5.3 %
Skew 1.3 stator slot pitch H= 7.0 %
6. CONCLUSION