tions in an electric machine may cause undesirable 'ripple' forces leading to acoustic noise. It is therefore important that the ripple fields and forces be computed ...
IEEE TRANSACTIONS ON MAGNETICS. VOL. 3 1 . NO. 3, MAY 1995
I472
Accurate Computation of 'Ripple Solutions' on Moving Finite Element Meshes Igor Tsukerman Department of Electrical and Computer Engineering, University of Toronto, Toronto, Ontario, M5S 1A4, Canada Abstract-The 'ripple field', i.e. the variation of the magnetic field for different positions of a moving part of an electric machine, can be computed accurately using the 'sliding surface' mesh model. The spurious ripple solution for this model is by a factor of O(h'/2) lower than for the models involving mesh distortion (h is the mesh size). Theoretical error analysis and numerical results are presented.
I. INTRODUCTION Variation of the magnetic field for different rotor positions in an electric machine may cause undesirable 'ripple' forces leading to acoustic noise. It is therefore important that the ripple fields and forces be computed accurately. The accuracy of Finite Element (FE) calculations of ripple solutions is related to the chosen method of movement modelling (see section I11 and the survey [4]). The objective of this paper is to compare the accuracy of ripple solutions for different types of moving FE meshes. Theoretical error estimates will be obtained in the form E = O(h"), meaning that c1ha 5 6 5 c2ha for some constants c1, c2 independent of the FE mesh size h. Although numerical values of these constants are not usually known, estimates of this kind do indicate how quickly the numerical error decreases when the FE mesh is refined. 11. FORMULATION
For purposes of theoretical analysis we shall consider a linear magnetostatic problem. However, numerical examples with nonlinearity and eddy currents will be considered in section V. The magnetostatic equation for the one-component magnetic vector potential A in a 2D domain with the reluctivity v(e,y, V A ) and the current density J ( z l y) is
divvVA
=
-J
(1)
with appropriate boundary conditions. For a rotating electric machine, we are interested in computing the 'ripple potential', i.e. the difference SA between the vector potentials A I , A2 for two nearby rotor positions (Fig. 1):
6A = Ai
-
A2
(2)
Fig. 1. Flux lines in a synchronous motor for two nearby rotor positions. Open circuit regime.
be emphasized that 6 A is meant to be a finite difference, not an infinitesimal variation/differential. When the equation (1) is solved using FEM on a mesh with a characteristic discretization size h, one obtains approximate numerical solutions Alh, A a h for A I , Az. The numerical value of the ripple potential (2) is bAh
=
Alh
-
A2h
(3)
The objective is to estimate the accuracy of the numerical ripple potential bAh, i.e. to compare 6Ah (3) with SA (2). One of the applications of this analysis is to the computation of radial ripple forces. 111. MODELLING THE MOTION
A. Methods The difficulty of modelling the motion lies in matching the solutions in the rotor and stator subdomains. The known techniques can be divided in two major groups. Hybrid methods are based on a combination of FEM with Fourier series [l],[3] or integral equations [9] in the air gap, while 'pure FE' models employ FE meshes in the whole domain, including the air gap [a], [7], [8], [12]. We shall consider only pure FE models and, even more specifically, the Sliding Surface (SS) [7] and the moving band [8] type methods with first order triangular elements.
These two positions are assumed to be close enough, so that SA is much smaller than A I , A2. However, it should
B. Moving Band Model
Manuscript received July 7, 1994. The work was supported in part by the Natural Sciences and Engineering Research Council, Canada, and by GE Canada,Inc., Peterborough, Ontario, Canada.
A 'moving band' [8] is a layer of air gap elements that are distorted during the motion. One side of this layer belongs to the moving part, and the other side is stationary. When
0018-Y464/Y5$04.00GI 199.5 IEEE
1473
ROTOR Fig. 2. Sliding surface in the air gap.
the rotor moves, the shape of the elements in the band changes, and at certain moments of time nodes on the sides of the band are re-connected to avoid significant distortion of the elements. Different variations of this approach are possible: the distortion can be spread over a few layers of elements or even over the whole air gap which should be re-meshed from time to time as suggested in [2]. Consider a hypothetical case of a smooth (toothless) axisymmetric stator with no stator load. Clearly, because of the axial symmetry, there is no ripple field associated with rotor motion, i.e. for any two rotor positions A1 = A;! E A , and 6A E 0. However, due to the well known estimate of the FE analysis [5], [6], the norm of the numerical ripple solution for the first order elements will be
II 6Ah IIE
= O ( h ) II A
IIW;(n)
where W; is the Sobolev space [5], and
(sa
-
11 .
Fig. 3. 'Layer Function' in a Fragment of the Air Gap.
define separate boundary value problems in the rotor and stator subdomaind. Consider the rotor subdomain for definiteness. In order to obtain a FEGalerkin formulation with the homogeneous Dirichlet condition on SS, we decompose the FE functional space Wh spanned by the basis functions in the rotor in two subspaces:
W h =
wp
@
W p
where the subspaces Wf) and Wfn) correspond, respectively, to the nodes on SS and to the inner nodes of the rotor subdomain. The vector Ah is split accordingly: (7)
(4)
is the energy
112
v V A . V A d S ) . Thus the numernorm 11 A I I E = ical ripple solution in this case is completely spurious and its order of magnitude is O ( h ) .
(6)
A t " ) is spanned by the basis functions for inner nodes in the rotor, and A t ) is a "layer function" (Fig. 3) spanned by the FE basis functions ?,bk for the nodes on SS:
Ai) =
AI(C)?,bk
(8)
node k E SS
C. Sliding Surface Model
A Sliding Surface (SS) [7] with equidistant nodes is generated inside the air gap (Fig. 2). As the pairs of coinciding nodes on SS should be matched, only discrete rotor shifts (by an integer number of the SS subdivisions) are allowed. However, this is not a significant restriction because, first, the discretization of SS is normally sufficiently fine, and, second, changes in the solution corresponding to movements much smaller than the mesh size cannot be obtained accurately anyway. On the other hand, the SS model has an important advantage of avoiding mesh distortion during the motion. As shown in section IV, this is essential for accurate computation of ripple quantities.
Iv. ERRORANALYSIS FOR SLIDING SURFACE MODEL Let A s s ( 4 ) be the ( a priori unknown) distribution of the vector potential on SS. Then the equation (1) with the boundary condition
Note that functions marked with a superscript "(in)" are zero on SS, while functions with the superscript "(1)" are nonzero only in one layer of elements adjacent to SS. With this notation, the FEGalerkin form of the problem with the homogeneous Dirichlet boundary condition is
( v V A ! " ) , VA',)
+ ( v V A t ) ,VA'h)
= -(J, A',) (9)
where (,) is the conventional inner product in Lz(f2);A i is an arbitrary trial function in Wp). This is in fact the standard weak formulation which is (often implicitly) used for inhomogeneous Dirichlet conditions. Due to the assumed linearity of the problem, equations (l),(9) are also valid for variations of the exact and numerical vector potentials, respectively. Thus (for SJ = 0)
diu v V ( 6 A ) = 0 (vVSA(h'n),VA',) = - ( v V 6 A t ) , VA',)
(10)
(11)
'For the rotor subdomain,the problemis consideredinthe reference frame moving with the rotor.
1474
For the 'ripple potential', the boundary condition on SS is
SA = SAss(4)
(12)
the maximum-norm error for Ah can be expressed via the maximum-norm approximation error and the solution error in a weaker norm [lo], e.g.:
differ from the exact values, so
SA(,') = SA(') +
+ II A - Ah llL2(nl)) E(')
(13)
Therefore (11) is rewritten as
=
( v V SA(hin),VA',)
- ( v V SA('),VA',)
- (vVc('), VA;)
(14)
SA!") can be split in two parts corresponding to the two respective terms in the right hand side of (14):
(22)
Both norms in the right hand side of (22) are known to be on the order of O(h2) [5], [6], so that
11 A
-
VA' E W Q l )
-
1
Ah IIL,(%)
5 O(h2)log5; 11 A
I\W&(%)
and therefore
II € ( I )
IILa(SS)
I I1 ( A 2 - A(',1
E
II SAP -
IlLa(SS)
SA(') IIL2(SS) 5
+ II ( A 2 - 4)) IlL2(SS) L 1
I W 2log5; ) II A Ilw&(nl)
where
( v V S A ~ )VA',) ,
= -(vVSA('), VAL)
(vVSAY), V A L ) = - ( v V d ' ) , VA',)
(16) (17)
+
Note that the sum SAf) SA(') is the FE solution corresponding to the boundary value problem (lo), (12) for the theoretical variation SA, while SA?) is due to the fact that the node values of the numerical solution on SS are not exact. The standard estimate of the FE analysis [5], [6] shows that
11 SAT) + 6A(')- SA 1 1 ~5
ch
11 SA Ilw;(n)
(18)
It is now left to estimate 6 A t ) . Taking A', = SA?) in (17) and using the Schwartz inequality, we obtain
(vVSA~V ) ,S A t ) ) = -(vVc('), V S A Y ) ) 5
(23)
(24)
Collecting all intermediate expressions (7), (13), (15), (18), (20), (21), (23), we obtain the final result:
11 bAh - SA I I E = 11 6At") + SA!) 5 11 SA?)
+ 11 e(') IIE
- SA 1
1 ~5
+ &A(')- 6 A 1 1 +~ 11 &A(,')1 1 +~ 5 ch 11 SA IlWi(s2) + 2 11 e(') llE
I
1
I O ( h ) II SA Ilw,a(n) + O ( h 3 / 2 1% ) 5; I1 A Ilw.&(nl) (25) The last term in (25) can be viewed as a 'spurious ripple solution' because it is proportional to the norm of the full solution A as opposed to the norm of the vuriution SA. However, this 'spurious' term has O(h1l2)times lower order of magnitude with respect to the mesh size than a similar spurious term (4) for models involving mesh distortion.
V. NUMERICAL EXAMPLES
5 11 6")
IlE
. 11 S A t ) IIE
(19)
3 11 b A t ) IIE 5 1) f(') llE (20) The norm 11 c(') in (20) can be estimated via the trace of d') on SS. Since d') by definition is nonzero only in the layer of thickness O(h) adjacent to SS, the gradient of E(') inside the layer is O ( h - l ) €('I, and
a 11 6'')
I)E
= O(h-1'2)
11 6")
I)La(SS)
(21)
Schatz and Wahlbin [lo] obtained the following "interior maximum norm estimate" for the error of the FE solution. Let R1 be a subdomain of the whole computational domain R and Rz be a subdomain of s21, i.e. R2 CC 01 CC Rz. (For our problem, RI can be the air gap or its part.) Then
As a model problem, we consider computation of the Mdial ripple force acting on the stator of a 20 pole, 2250 kVA synchronous motor [13], [14]. The FE mesh for a single pole had 2856 nodes and 5186 elements with 199 nodes on SS. Magnetostatic and time dependent eddy current problems (eddy currents in the field and amortisseur windings) for the open circuit regime were considered. The first test is for the hypothetical case of a smooth stator (no teeth and slots) when theoretically the ripple force is nil. While the numerical ripple force in the SS model is virtually zero, there is a spurious ripple force for the method with moving band distortion. If the elements of the moving band are allowed to get shifted by about three mesh sizes without re-meshing (Fig. 4), the spurious ripple solution reaches 1% of the full radial force (which is an appreciable part of the ripple force). If the distortion of the moving band is kept within a limit of one mesh size,
1475
.4-
2.
Relative
RippkForce,% 0 .
-2
0
-2.
2
Relative Rotor Displacsmmt
Fig. 4. Spurious Ripple Force vs. Relative Rotor Displacement for the Moving Band Method.
-0
0.001
0.002
0.003
0.004
Time [seconds] Fig. 6. Ripple Force vs. Time for Sliding Surface Model. Eddy TC model, 0 SC model. Current Case.
-
0.5.
Relative Ripple Force, 5%
developed jointly by the University of Toronto (E. Sidoriak, I. Tsukerman and J.D. Lavers) and GE Canada, Inc. (K. Weeber, K. Nguyen, H. Karmaker). I am especially grateful to Evan Sidoriak for his help and advice on software issues.
O'
-0.5
-0
0.002 Time [seconds]
0.001
0.003
Fig. 5 . Ripple Force vs. Time for Sliding Surface Model. Magneto- TC model, 0 SC model. static Case.
then the spurious ripple force remains small (about 0.1% of the total force). However, one should be aware of the spurious ripple force effect when applying methods with mesh distortion and re-meshing. The second test is for the real machine geometry (Fig. 1). Two SS mesh models have been generated: tooth centered (TC) for the rotor d-axis positioned against the middle of a stator tooth, and slot centered (SC) for the d-axis against a slot. The steady state value of the ripple force should be the same for the T C and SC models; this is a consistency check for the simulations [13], [14]. Numerical results for TC and SC starting positions are indeed very close in both magnetostatic (Fig. 5) and eddy current (Fig. 6 ) problems. (In these particular simulations, the difference between the force waveforms of Fig. 5 and 6 is caused by the amortisseur bar currents.) More discussion of the magnetostatic case can be found in [13]. VI. CONCLUSIONS Error analysis of the numerical 'ripple solution' on a moving FE mesh indicates a 'spurious' component which is related to movement or distortion of the mesh. For the sliding surface mesh model this spurious component is shown to be substantially smaller (by the order of O ( h l / ' ) )than for models involving mesh distortion. ACKNOWLEDGMENTS Software for simulation of time-dependent fields and eddy currents in synchronous machines (section V) was
REFERENCES [l] A.A. Abdel-Razek, J.L. Coulomb, M. Fe'liachi, J.C. Sabon-
nadibre, "The calculation of electromagnetic torque in saturated electric machines within combined numerical and analytical solutions of the field equations," IEEE Trans. Magn., 17(6),pp.32503252,1981. [2] S.R.H. Hoole, "Rotor motion in the dynamic finite element analysis of rotating electricalmachinery," IEEE Trans. Magn., 21(6), pp.471474, 1985. [3] Ki-sik Lee, M.J. DeBortoli, M.J. Lee, S.J. Salon, "Coupling finite elements and analytical solution in the airaap of electric machines," IEEE Trans. Magn., 2 7 ( 5 ) , pp.3955-g957,1991. [4] J.D. Lavers, "Electromagnetic field computation in power engineering," IEEE Trans. Magn., 29(6), pp.2347-2352,1993. [5] J. T. Marti, Introduction to Sobolev Spaces and Finite Element Solution of Elliptic Boundary Value Problems, Academic Press, 1986. [6] L.A. Oganesian, and L.A. Rukhovets, Variatsionno - Raznostnye Metody Resheniia Ellipticheskikh Uravnenij (VariationalDifference Methods for the Solution of Elliptic Equations): Izd-vo AN Armianskoi SSR,Yerevan, USSR [in Russian]. [7]T.W.Preston, A.B.J.Reece, P.S.Sangha, "Induction motor analysis by time-stepping techniques," IEEE Trans. Magn., 24(1), pp.471474, 1988. [8] N. Sadowski, Y.Lefevre, M. Lajoie-Mazenc, J. Cros, "Finite element torque calculation in electrical machines while considering the movement," IEEE Trans. Magn., 28(2), pp.1410-1413,1992. [9] S.J. Salon, J.M. Schneider, "A hybrid finite element - boundary integral formulation of the eddy current problem," IEEE Trans. Magn., 18(2), pp.461466,1982. [lo] A.H. Schatz and L.B. Wahlbm, "Interior maximum norm e% timates for finite element methods," Math. of Comp., 31(138), pp.414442,1977. [I11 R. Scott, "Optimal L , estimates for the finite element method on irregular meshes," Math. of Comp., 30, pp.681-697,1976. [12] I.A. Tsukerman, "Overlapping finite elements for problems with movement," IEEE Trans. Magn., 28(5), 2247-2249,1992. [13] I. Tsukerman, J.D. Lavers, A. K o ~ r a d "Using , complementary formulations for accurate computations of magnetostatic fields and forces in a synchronous motor," IEEE Trans. Magn., 30(5), pp.3479-3482,1994. [14] I. Tsukerman, J.D. Lavers, A. Konrad,K. Weeber, H. Karmaker, "Finite element analysis of static and time-dependent fields and forces in a synchronous motor," Proceedings of the International Conference on Electrical Machines, Paris,1994, 2, pp.27-32.
