- Email: [email protected]
A finite element solution for periodic eddy current problems in hysteretic media
Journal of Magnetism and Magnetic Materials 160 (1996) 96-97
~ Journalof magnetism J H and magnetic materials
,~
ELSEVIER
A finite element solution for periodic eddy current problems in hysteretic media O. Bottauscio a, M. Chiampi b,,, M. Repetto b a lstituto Elettrotecnico Nazionale Galileo Ferraris, Torino, Italy b Dip. lngegneria Elettrica lndustriale, Politecnico di Torino, c.so Duca degli Abruzzi 24, 10129 Torino, Italy Abstract
A finite element formulation for hysteretic media, under periodic supply, is presented. The computational scheme, including eddy currents, is developed in the harmonic domain and uses the fixed point technique and scalar Preisach model in reversed fashion. The method is validated by comparison with experimental results. Keywords: Eddy currents; Finite element method; Hysteresis model 1. I n t r o d u c t i o n
This paper presents a finite element formulation for the solution of time periodic problems including hysteresis and eddy currents. Under the 2D hypothesis, the field problem is formulated in terms of the magnetic vector potential and the nonlinearity is tackled by the fixed point (FP) technique. The hysteretic behaviour is modelled using the classical scalar Preisach scheme [1] in the reversed fashion, since the field problem provides the waveform of the flux density B ( t ) . The periodic problem is solved in the frequency domain; then the harmonic solution is translated into the time domain to handle the nonlinearity; finally, a fast Fourier transform (FFT) is used to update field equations for the next iteration. The procedure is applied to the analysis of magnetic sheets, to predict the dynamical hysteresis loops. The results are in good agreement with experiment. 2. F i n i t e e l e m e n t s o l u t i o n s c h e m e a n d n u m e r i c a l r e s u l t s
The 2D eddy current problem in a magnetic sheet is governed by
VX[~'(VXA)]=J-o-~-
OA
~ r OA
t +--I.--dO s .~ at
(1)
together with the stated boundary conditions. In Eq. (1), J is the imposed current density, A the magnetic vector potential, o" the conductivity, I2 the cross section of the magnetic material with area S and H = ( ( B ) is the hysteretic. The nonlinearity is handled by the FP technique
* Corresponding author. Email: [email protected]; fax: + 39-11-5647199.
which splits ~ into a linear term and a residual R ( ( ( B ) = UTB + R). Coefficient uv depends on curve ff and R has to be evaluated iteratively. The FP method is attractive because it ensures convergence, does not require a smooth curve and can be used in conjunction with polydrome characteristics. Under periodic conditions, Eq. (1) is handled in the harmonic domain, introducing a set of linear equations: _ -- j wno-A~(,~) vx (~Tv ~_k l = ~'~ --
X
--(n)~
.
•
o"
+jwn--f
SJa-
,,,)
A k d/2-
V×
R(,, ) _k
],
(2) where o~ is the fundamental angular frequency, n the (n) harmonic order, k the iteration number and A k is a phasor. A weak form of Eq. (2) is solved by the finite element (FE) method using linear shape functions. Since the FP scheme only changes the right-hand side (RHS) of Eq. (2), the stiffness matrix is only computed once. At each iteration, the field equations provide the harmonic components of nodal values of A and the magnetic flux density on each element. The time behaviour of B is obtained by an inverse FFT. To compute residual nonlinearity R k, H waveforms have to be deduced from B evolution using the Preisach model in the reversed fashion. Knowing the time evolution of R k, the FFT provides the harmonic distribution of the residual nonlinearity and allows updating of the RHS of Eq. (2). The procedure, sketched in Fig. 1, is iterated until convergence. The heart of the computational scheme lies in the use of the Preisach model in the reversed fashion. Under periodic supply, the wipe out property of the Preisach model can be exploited in order to determine a starting point from which the next evolution can be predicted. In
0304-8853/96/$15.00 Copyright © 1996 Elsevier Science B.V. All rights reserved. PI1 S0304-8853(96)00124-2
O. Bottauscio et al. / Journal of Magnetism and Magnetic Materials 160 (1996) 96-97 [ Spatial discretization~
I
--.---v~
I
[=-0
-~o--
['=3T/8
•
97
t=- T I 8
--~---
['=- TI4
[=-TI2
1,51
field solutions ]
=AAAA&&&A&&&m,&&&A&A&A&&&&A&&&&&&&&A&AA&'
A(o~) B(co) I inverseFFY ]
~'
!
1,0
B(t) l Inverse Preisach ]
0,5
0,0
No
"w e..
.-0,5 -1,0
-1,5
Fig. 1. Flowchart of the procedure.
0,0
/ \
............
~
L
0,2
0,4
,,' ',,.
L 0,6
h 0,8
1,0
Laminationthickness [mm]
fact, the instant when B reaches its peak value, univocally defines the staircase interface in the Preisach plane and the corresponding ( H , B ) point lies on the normal magnetization curve. In successive instants, an iterative search for the H value is employed, starting from the known configuration of the Preisach triangle at the previous step. A secant~bisection procedure is particularly efficient for this purpose [2]. Once H is obtained at the actual time, the corresponding staircase interface becomes the starting point for the next evolutions. Since the Preisach scheme is used in the reversed fashion, extension to the moving Preisach model requires only minor modifications. The magnetization is computed by numerical integration of the distribution function over the Preisach triangle. This operation, included in the iterative process, is repeated several times and so it is essential to reduce its computational burden. Unfortunately, the distribution function often shows very high and narrow peaks which are not well handled by the usual multidimensional quadrature rules; on the other hand, a too refined numerical integration over the whole Preisach plane leads to unacceptable computational costs. An adaptive quadrature rule based on the Kronrod [3] scheme has been found to be the best compromise between accuracy and computational burden. The procedure is accelerated by preliminary integration of the distribution function over each element of a nonuniform mesh created on the Preisach triangle.
0,5 '~
0,0
=z~ -0,5
Fig. 3. Instantaneous magnetic flux density distribution. The distribution function q~ is expressed in the form ¢(c~,/3)= O ( a ) 0 ( - / 3 ) . Thus, linking N experimental variations of the magnetization over one branch of the main hysteresis loop with the integral of ¢ over the Preisach plane, a bilinear system of equations is obtained [4], where the unknowns are N sampled values of function tp. Function ~ is then interpolated linearly in order to compute the distribution function. The proposed method has been applied to the analysis of a 1 mm thick FeSi 2 wt% lamination, where a sinusoidal magnetic flux is imposed. Because of saturation, the imposed sinusoidal flux creates a harmonic spectrum of local magnetic quantities. The amplitude of this spectrum is unknown a priori and was estimated by a trial and error procedure. In all tested cases, a 32 harmonics spectrum has been shown to be sufficiently accurate. The computation was developed considering different values of the imposed flux and supply frequency (from 5 to 1500 Hz). A comparison with experiments was performed on the dynamic hysteresis loops. Computed and experimental results are in good agreement, as shown in Fig. 2, where a 50 Hz sinusoidal flux which gives rise to a mean magnetic flux density of 1.5 T is considered. The corresponding instantaneous distribution of the magnetic flux density over the lamination thickness is presented in Fig. 3, where five different instants in the semiperiod T / 2 are examined. It is worth noting that, at t = T/4, due to saturation, the skin effect is negligible. The figure also shows a significant phase difference in B waveforms at different points along the lamination thickness, which causes an increase of the magnetic losses.
References
-1.0 -1.5 -750''-5~''
'-.2'50 ' " " 0 . . . .
2'50. . . .
500 . . . .
750
H [A/m]
Fig. 2. Measured and computed hysteresis loop.
[l] I.D. Mayergoyz, Mathematical Models of Hysteresis (Springer, New York, 1991). [2] O. Bottauscio et al., INTERMAG '95, San Antonio, Texas. [3] R. Piessens et al., Quadpack (Springer, Berlin, 1983). [4] G. Kadar, J. Appl. Phys. 61 (1987) 4013.
~ Journalof magnetism J H and magnetic materials
,~
ELSEVIER
A finite element solution for periodic eddy current problems in hysteretic media O. Bottauscio a, M. Chiampi b,,, M. Repetto b a lstituto Elettrotecnico Nazionale Galileo Ferraris, Torino, Italy b Dip. lngegneria Elettrica lndustriale, Politecnico di Torino, c.so Duca degli Abruzzi 24, 10129 Torino, Italy Abstract
A finite element formulation for hysteretic media, under periodic supply, is presented. The computational scheme, including eddy currents, is developed in the harmonic domain and uses the fixed point technique and scalar Preisach model in reversed fashion. The method is validated by comparison with experimental results. Keywords: Eddy currents; Finite element method; Hysteresis model 1. I n t r o d u c t i o n
This paper presents a finite element formulation for the solution of time periodic problems including hysteresis and eddy currents. Under the 2D hypothesis, the field problem is formulated in terms of the magnetic vector potential and the nonlinearity is tackled by the fixed point (FP) technique. The hysteretic behaviour is modelled using the classical scalar Preisach scheme [1] in the reversed fashion, since the field problem provides the waveform of the flux density B ( t ) . The periodic problem is solved in the frequency domain; then the harmonic solution is translated into the time domain to handle the nonlinearity; finally, a fast Fourier transform (FFT) is used to update field equations for the next iteration. The procedure is applied to the analysis of magnetic sheets, to predict the dynamical hysteresis loops. The results are in good agreement with experiment. 2. F i n i t e e l e m e n t s o l u t i o n s c h e m e a n d n u m e r i c a l r e s u l t s
The 2D eddy current problem in a magnetic sheet is governed by
VX[~'(VXA)]=J-o-~-
OA
~ r OA
t +--I.--dO s .~ at
(1)
together with the stated boundary conditions. In Eq. (1), J is the imposed current density, A the magnetic vector potential, o" the conductivity, I2 the cross section of the magnetic material with area S and H = ( ( B ) is the hysteretic. The nonlinearity is handled by the FP technique
* Corresponding author. Email: [email protected]; fax: + 39-11-5647199.
which splits ~ into a linear term and a residual R ( ( ( B ) = UTB + R). Coefficient uv depends on curve ff and R has to be evaluated iteratively. The FP method is attractive because it ensures convergence, does not require a smooth curve and can be used in conjunction with polydrome characteristics. Under periodic conditions, Eq. (1) is handled in the harmonic domain, introducing a set of linear equations: _ -- j wno-A~(,~) vx (~Tv ~_k l = ~'~ --
X
--(n)~
.
•
o"
+jwn--f
SJa-
,,,)
A k d/2-
V×
R(,, ) _k
],
(2) where o~ is the fundamental angular frequency, n the (n) harmonic order, k the iteration number and A k is a phasor. A weak form of Eq. (2) is solved by the finite element (FE) method using linear shape functions. Since the FP scheme only changes the right-hand side (RHS) of Eq. (2), the stiffness matrix is only computed once. At each iteration, the field equations provide the harmonic components of nodal values of A and the magnetic flux density on each element. The time behaviour of B is obtained by an inverse FFT. To compute residual nonlinearity R k, H waveforms have to be deduced from B evolution using the Preisach model in the reversed fashion. Knowing the time evolution of R k, the FFT provides the harmonic distribution of the residual nonlinearity and allows updating of the RHS of Eq. (2). The procedure, sketched in Fig. 1, is iterated until convergence. The heart of the computational scheme lies in the use of the Preisach model in the reversed fashion. Under periodic supply, the wipe out property of the Preisach model can be exploited in order to determine a starting point from which the next evolution can be predicted. In
0304-8853/96/$15.00 Copyright © 1996 Elsevier Science B.V. All rights reserved. PI1 S0304-8853(96)00124-2
O. Bottauscio et al. / Journal of Magnetism and Magnetic Materials 160 (1996) 96-97 [ Spatial discretization~
I
--.---v~
I
[=-0
-~o--
['=3T/8
•
97
t=- T I 8
--~---
['=- TI4
[=-TI2
1,51
field solutions ]
=AAAA&&&A&&&m,&&&A&A&A&&&&A&&&&&&&&A&AA&'
A(o~) B(co) I inverseFFY ]
~'
!
1,0
B(t) l Inverse Preisach ]
0,5
0,0
No
"w e..
.-0,5 -1,0
-1,5
Fig. 1. Flowchart of the procedure.
0,0
/ \
............
~
L
0,2
0,4
,,' ',,.
L 0,6
h 0,8
1,0
Laminationthickness [mm]
fact, the instant when B reaches its peak value, univocally defines the staircase interface in the Preisach plane and the corresponding ( H , B ) point lies on the normal magnetization curve. In successive instants, an iterative search for the H value is employed, starting from the known configuration of the Preisach triangle at the previous step. A secant~bisection procedure is particularly efficient for this purpose [2]. Once H is obtained at the actual time, the corresponding staircase interface becomes the starting point for the next evolutions. Since the Preisach scheme is used in the reversed fashion, extension to the moving Preisach model requires only minor modifications. The magnetization is computed by numerical integration of the distribution function over the Preisach triangle. This operation, included in the iterative process, is repeated several times and so it is essential to reduce its computational burden. Unfortunately, the distribution function often shows very high and narrow peaks which are not well handled by the usual multidimensional quadrature rules; on the other hand, a too refined numerical integration over the whole Preisach plane leads to unacceptable computational costs. An adaptive quadrature rule based on the Kronrod [3] scheme has been found to be the best compromise between accuracy and computational burden. The procedure is accelerated by preliminary integration of the distribution function over each element of a nonuniform mesh created on the Preisach triangle.
0,5 '~
0,0
=z~ -0,5
Fig. 3. Instantaneous magnetic flux density distribution. The distribution function q~ is expressed in the form ¢(c~,/3)= O ( a ) 0 ( - / 3 ) . Thus, linking N experimental variations of the magnetization over one branch of the main hysteresis loop with the integral of ¢ over the Preisach plane, a bilinear system of equations is obtained [4], where the unknowns are N sampled values of function tp. Function ~ is then interpolated linearly in order to compute the distribution function. The proposed method has been applied to the analysis of a 1 mm thick FeSi 2 wt% lamination, where a sinusoidal magnetic flux is imposed. Because of saturation, the imposed sinusoidal flux creates a harmonic spectrum of local magnetic quantities. The amplitude of this spectrum is unknown a priori and was estimated by a trial and error procedure. In all tested cases, a 32 harmonics spectrum has been shown to be sufficiently accurate. The computation was developed considering different values of the imposed flux and supply frequency (from 5 to 1500 Hz). A comparison with experiments was performed on the dynamic hysteresis loops. Computed and experimental results are in good agreement, as shown in Fig. 2, where a 50 Hz sinusoidal flux which gives rise to a mean magnetic flux density of 1.5 T is considered. The corresponding instantaneous distribution of the magnetic flux density over the lamination thickness is presented in Fig. 3, where five different instants in the semiperiod T / 2 are examined. It is worth noting that, at t = T/4, due to saturation, the skin effect is negligible. The figure also shows a significant phase difference in B waveforms at different points along the lamination thickness, which causes an increase of the magnetic losses.
References
-1.0 -1.5 -750''-5~''
'-.2'50 ' " " 0 . . . .
2'50. . . .
500 . . . .
750
H [A/m]
Fig. 2. Measured and computed hysteresis loop.
[l] I.D. Mayergoyz, Mathematical Models of Hysteresis (Springer, New York, 1991). [2] O. Bottauscio et al., INTERMAG '95, San Antonio, Texas. [3] R. Piessens et al., Quadpack (Springer, Berlin, 1983). [4] G. Kadar, J. Appl. Phys. 61 (1987) 4013.