Direct Consideration of Eddy Current Losses in Laminated Magnetic Cores in Finite Element Method (FEM) Calculations Using the Laplace Transform

: The following article presents a computation procedure that enables us to simulate the dynamic states of electric machines with a laminated magnetic core, with direct consideration of the eddy current losses. The presented approach enables a signiﬁcant reduction of the simulation process computational complexity. The veriﬁcation of the obtained data correctness is based on a detailed balance of energy and power in the investigated system. The correctness of the obtained results was also conﬁrmed by comparing them with the results included in norms that describe the losses in laminated sheets. The presented approach is based on expressing the equivalent permeability of transformer metal sheets by using RC or RL circuits. The impedances of these circuits are treated as the transmittance of Inﬁnite Impulse Response ﬁlters (IIR) of the Laplace s variable. In this form they are implemented in direct calculations of the dynamics of electric machines based on ﬁeld-circuital models, using the Finite Element Method (FEM). In this way, we present the method of including eddy current losses in laminated metal circuits of chokes or transformers, during calculations using the ﬁnite element method, with the IIR ﬁlter in the domain of the variable s of the Laplace transform. Eddy current losses are directly included in the calculation process. Therefore, they have a direct impact on the transient state waveforms. However, the use of the Laplace variable s caused an excessive increase in the number of state variables, and the overall computational efﬁciency of the presented method is sufﬁciently low so as to be used in the simulation process of electrical machine dynamic states with a relatively large number of elements in the FE Model.


Introduction
The prediction of losses in the magnetic cores of electric machines is an important issue in their design. At present, these machines are more and more frequently powered by inverters [1,2]. The non-sinusoidal shape of the voltage supplying electric machines significantly increases losses in their laminated cores, which makes this issue even more important.
These issues are of particular importance in the case of designing critical drive systems. Proper prediction of power losses in a machine core enables, even at an early design phase, protection against the results of excessive heating of the system from the results of thermal energy losses in the core. Critical drive systems are designed in such a way that they can continue operation also in emergency states. This may potentially lead to permanent excesses of the nominal values of machine currents, the significant influence of which, on the power of losses in the core, must be taken into consideration even at the design phase. W Fe = k h · B α + k e · f 2 · B 2 + k a · f 1.5 · B 1.5 (1) where k h is the coefficient referring to losses from eddy currents, k e is the coefficient referring to losses from eddy currents, and k a is the coefficient referring to overload losses (anomalous). The values and interpretation of these coefficients are described in detail in [16,17]. Typically, the coefficient α is equal to 2. Eddy current and excess losses can be grouped together and called eddy current losses: where f is the frequency. The Formula (2) works well with adequately chosen kh and ke coefficients, especially in the range of higher magnetic flux density ( 1 T). The following formula can be used to identify the coefficients k h and k e : Depending on the frequency of the sine waves of the magnetic flux density B, from Formula (3) different values of k h and k e coefficients are obtained subjected to the magnetic flux density. In many cases, electric machines are supplied with non-sinusoidal voltages from transducers. The losses in iron are increasing. However, their determination requires prior information on the shape of the applied voltages, as they result in an appropriate magnetic flux density B(t). In the works [11,12], it was justified that the maximum (peak) value of magnetic flux density, decisive for hysteresis losses, is determined by the average value of the supply voltage V avr . As the eddy current losses are dependent on the square of the magnetic flux density, the derivative d dt B(t), they can also be dependent on the effective supply voltage V RMS . Therefore, for any shape of the supply voltage, the iron losses of the machine can be determined from the formula: where η = V avr V 1,avr ; ε = V RMS V 1,avr are the relative quantities concerning the first harmonic. Other methods for determining the iron losses with distorted supply voltages are included in the works [10,[18][19][20]. The commercial Ansoft Maxwell software uses the Equivalent Elliptical Loop (EEL) model to directly calculate the hysteresis losses P hy : where: B m is the maximum value of magnetic flux density, and B is the current value of the magnetic flux density. Therefore, Formula (6) requires several calculation periods to determine the maximum value of B m . This process is described in [15]. By using the above-mentioned methods, it is possible to determine the value of the equivalent resistance R Fe , the inclusion of which in the electric machine equivalent scheme will enable us to consider the influence of the presence of losses during simulations of dynamic states. The R Fe should be included in the equivalent scheme in parallel to voltages that are derived from coupled fluxes for each phase of the machine winding. This was presented with the example of the SRM motor [21,22]. Another approach to consider the losses during dynamic simulation, by modifying the equivalent scheme, is presented in [23]. The parameters of an equivalent scheme of a transformer are adjusted based on the frequency response.
It should be emphasized that most of the methods of eddy current loss calculations used in the literature are inaccurate. In order to determine the losses in laminated magnetic systems, the equivalent conductivity of sheets is usually determined, followed by calculations using the finite element method [14,24]. However, this requires the assumption of a priori waveforms to calculate the equivalent conductivity of sheets. The inaccuracy is caused by the fact that the assumed spatial distribution of magnetic flux density does not correspond to the real one.
A very accurate method of considering the eddy currents in laminated sheets of magnetic cores was presented in [25,26]. This method considers the effects of eddy currents in the core without the need of separate modeling of all sheets. The solution of the Helmholtz (8) equations that describes the phenomena that take place in laminated sheets by means of series is used there. The component elements of the series are products of the assumed base (shape) functions dependent on the space coordinate perpendicular to the surface of sheets and the searched function dependent on time. These dependent on time functions occur directly in the assumed distribution of magnetic induction in a series and as derivatives of time in distribution into a series of magnetic field intensities.
After application of the finite element method in the plane parallel to the surface of sheets, the components of magnetic induction, dependent on time, are unknown elements. They are expressed by means of vector magnetic potentials and the potentials are calculated in an iterative process. Depending on the need to consider the occurrence of higher harmonic elements of harmonic waveforms of time, the number of the components of distribution into a series should increase. This inflates the number of unknown elements. Such distributions into series for induction and magnetic field strength should be created for each finite element. The equations are solved by means of the iterative Newton-Raphson method. Due to a very high number of unknown elements and the use of iterative processes, the method is very demanding in terms of computation time. In the case of complex structures that occur in the calculated electric machines, this requires a high number of finite elements. These computation difficulties may limit the use of this method.
The considerations presented in this article will serve to account for eddy current losses in laminated sheets of the magnetic core of electrical machines, directly during the calculations of dynamics, using Infinite Impulse Response filters (IIR) filters.
Thanks to a significant reduction of the computational complexity of the simulation process, in relation to the methods described in [25,26], the approach presented in this article enables the modeling of electric machines by using a high number of finite elements.
There are methods, presented, for example, in [27], that allow the consideration of the skin effect and the proximity effect in a quad of electric cables by means of the homogenization (unification) method. Thanks to that, consideration of the phenomena that take place separately in particular windings is avoided. Finite Element Analysis (FEA) in the domain of time with a homogenous winding requires the introduction of additional unknown elements and equations. This extension in the domain of time is analogous to induction impedance approximation dependent on frequency, by means of a ladder system that consists of constant resistances and inductances.
The article [28] makes use of the Lagendre polynomial for the purpose of describing flux induction distribution in a magnetic plate. On this basis, an equivalent Cauer scheme is derived for the purpose of modeling eddy current losses in a laminated sheet. A similar idea of deriving an equivalent scheme for the magnetic permeability of laminated sheets in the form of an l, R scheme and its application in field analysis by means of Finite Element Method (FEM) is a primary idea of this article.

Direct Synthesis of Equivalent Magnetic Permeability of Laminated Sheets in the Form of a Transfer Function of an IIR Filter
In the presented method, the complex permeability,μ, of a laminated magnetic core is derived. It can be expressed by impedance of the R, C system presented in Figure 1. In a similar manner, the complex reluctivity of a laminated core,ν = 1 µ , may be expressed through impedance of the R, L circuit presented in Figure 2. The imaginary part of the magnetic permeability corresponds to active power, the real part to reactive power. These circuits should be treated as IIR filters, the transmittance of which comprises the representation of the equivalent magnetic permeability of a laminated core in the domain of the complex variable s of the Laplace transform. The equivalent magnetic permeability can be derived from Figure 3 and Formula (8) [29,30]. In order to include the saturation of the magnetic circuit in these calculations, an iterative process was developed. A scalar magnetic hysteresis model according to Juhani Tellinen has been implemented to also include hysteresis losses [31].
In order to consider the magnetic field in the transformer metal sheet, it is necessary to solve Maxwell's equations in the coordinate system, as shown in Figure 3 [4,32]. We assume that the height of the plate in the x direction and its length along the y axis are infinitely large, and the magnetic flux density has only the y component. The magnetic field evokes a current that has a component in the x-axis direction and that is unevenly distributed in the z-axis direction with the density of I(z). Maxwell's equations are used then.
where µ is the magnetic permeability of iron plate and γ is the electric conductivity. According to the considerations in [29], the formula was obtained for the mean value of the magnetic flux density in the plate B avr and for the complex magnetic permeabilityμ: where p = j · ω · µ · γ, while γ is the material conductivity and µ is the magnetic permeability.
The Equations (9) and (10) are valid for mono-harmonic waveforms of the given angular frequency ω. In order to obtain equations for any transient waveforms in time, it is necessary to substitute j · ω with the Laplace variable s. This substitution enables a shift from 'details', i.e., mono-harmonic waveforms, to 'generality', i.e., transient waveforms. Therefore, it is not always easy to perform it. However, it is used, for example, in the theory of digital filters, in order to derive equations for the transmittance of Butterworth or Czybyszew filters. Then, it is necessary to assume p = √ s · µ · γ.
The magnetic permeability in Equation (10) assumes an operational form, in which there occurs √ s. In order to transform this permeability to the form of a rational function of Laplace's variable s, it is possible to make use of the methods connected with fractional derivatives. When deriving Equation (11), an inverse Laplace transform [29] was used.
According to [29], the equivalent magnetic permeability of sheets from the Formula (10) can be represented in the form of the relationship: The equivalent magnetic permeability can then be replaced with the impedance of the circuit in Figure 1. For the purpose of practical calculations, the infinite number of R,C series circuits connected in parallel is reduced to the number of n. The remaining part, beginning with n + 1, is replaced with one equivalent RC branch. The parameters of this branch are selected with optimization methods that provide the best approximation in the given angular frequency range. The calculations are done with complex numbers. The unknown elements are the investigated parameters R and C of the parallel equivalent branch. What is searched here is the minimum of the sum of absolute values from the differences in the magnetic permeability from Equation (10) and impedance of the system from Figure 1, limited to n of parallel branches for the selected angular frequencies from the selected range from several Hz to 5000 Hz. In such a case, the classical optimisation algorithms for a convex problems are sufficient. When considering a wider range of angular velocities, it seems that it would be reasonable to apply the newest algorithms presented in [33,34]. The articles show the problem of identification of guaranteed passive rational models for linear distribution structures. The identification is pursued in two separate steps: In the first step, the Vector Fitting algorithm is used for identification of the poles, then the identification of residues is formulated as a convex programming problem with constraints on the single terms of the pole-residue expansion. The equivalent magnetic reluctanceν = 1 µ is in the denominator of Formula (11). It can be presented in the form of a formula:ν and in the form of the R,L impedance from Figure 2 Sample non-linear iron magnetization characteristics were used for calculations in the formula: where a 1 = 398; a 3 = 30; a 9 = 55. The condition of using a given formula as a magnetic characteristic is the one-uniqueness. In Figure 4 the difference between the curves marked as 1 and 2 for the assumed value of the magnetic flux density B determines the hysteresis losses. It presents the implementation of Tellinen's hysteresis model in the time calculations [29,35].

Inclusion of Eddy Current Losses in the Laminated Magnetic Circuit in Transient Finite Element Method Calculations
The transient state of the electric waveforms taking place in the choke should be calculated. A winding is wound on the magnetic core, which is switched on to a sinusoidal voltage. The magnetic core is made of laminated sheets, in which heat is released due to the eddy current losses. The system cross-section through the plane (x-y) is considered; therefore the symmetry along the coordinate axis z is assumed (Figure 7).
It is therefore necessary to calculate the magnetic flux density, B(x, y), and the magnetic field strength, H(x, y). A detailed derivation of the equation that describes the analyzed problem is presented in [36].
The equation that needs to be solved by means of the Finite Element Method assumes the following form: ∂ ∂x where µ x and µ y are the magnetic permeabilities.
It should be emphasized that, in this case, one is dealing with a magnetic anisotropic environment, thus µ x = µ y . According to the evidence presented in [37], it is necessary to apply rectangular coordinates (x, y) connected with the environment. This manifests in the magnetic permeability µ x along the x coordinate system and the magnetic permeability µ y along the coordinate system of y. The current density J z should be dependent on the current flowing through the windings, which is a state variable. The magnetic permeability should include the eddy current losses in the laminated magnetic circuit, as described above. The cross-section x-y of the coil wound on the magnetic core is assumed, as shown in Figure 7 [31]. The coil is surrounded by air, and the whole system is contained in a rectangular area. We assume that the magnetic field does not escape outside the rectangle, therefore the zero Dirichlet Boundary Condition is imposed on the vector potential A z [38]. To apply the weak Galerkin formulation method, Equation (14) is multiplied by any variation of the vector potential σA, using the following formula: and then integrates over the entire section area [19,39]. When integrating the right side of the formula (15), the Gauss theorem is used: where V is any volume limited by the surface S and vector dS is the vector of normal to elementary surface dS, directed outside the area V. Zero Dirichlet conditions for the vector potential A were assumed at the boundary S. Therefore, the change of this potential, i.e., the variation at this boundary S, is also equal to zero. The first part of the right side of the Equation (15) is therefore zero, and the weak Galerkin formulation for Equation (14) can be noted as: where V is the x-y cross-section of the choke shown in Figure 7. Formula (14) was used here. The left side of this formula, with the assumed symmetry to the axis z, is equal to: div( 1 µ · grad( A)) = J z . The system of equations of the finite element method is created by adding up the contributions from each finite element. Local coordinates are used [5]. The equivalent of the contribution of the first part of Formula (17) [5,31] to the general equation of the finite element method derived from the triangular finite element with linear shape functions is given by: where A i , A j , A m are the searched magnetic potential values in triangle nodes and coefficients. The coefficients b i , b j , b m and c i , c j , c m depend on the difference between the global coordinate nodes of the Finite Element Method [5].
To determine the equivalent magnetic reluctance of the laminated core sheets,ν = 1 µ , the formula (12) is used. The magnetic reluctanceν(s) is obtained in the form of a rational function of the Laplace variable s, i.e., in the form of a transfer function of the IIR analogue filter. The system of equations obtained using the Finite Element Method is then multiplied by the IIR filter denominator to find the solution. Thanks to this multiplication, theν(s) filter of the variable s will apply globally to all finite elements. This operation causes most of the computational costs and numerical difficulties in the presented procedure.
Formula (18) contains the essence of the implementation of the IIR filter, representing the equivalent magnetic permeability of laminated sheets, µ(s), for the Finite Element Model. The second segment of Equation (17), containing the current density J z should also be included in these calculations. This is the current density evenly distributed in the coil cross-section, caused by the coil current I c . Therefore, the unknowns are: vector potentialū = A z and current I c . The additional voltage equation should be created for the coil circuit, which is given by the Equation where ρ w is the winding density in the cross-section S, L rs and R s are the leakage inductance and winding resistance, and E(t) is the winding supply voltage.
In the case of a larger number of windings, a correspondingly larger number of additional equations must be taken into account. The alternative approach presented at this point, using the Laplace variable s to represent the equivalent magnetic permeability of transformer sheets of the magnetic circuit, will consist in multiplying each of the Equation (18) by denominators of rational functionsν(s) in the Laplace variable domain. The procedure during the implementation of theν(s) filter for the Finite Element Method can, therefore, be called global. If different saturation levels of the magnetic circuit are included for different finite elements, then the presented procedure requires multiplication of the equations by the product of the IIR filter denominators, which can lead to further numerical difficulties.
The presented procedure causes that individual unknowns are to be multiplied by the appropriate degrees of the Laplace variable s, up to the coefficient s N f ilter −1 , where N f ilter is the degree of the filter.
The unknown creates a vector (the vector potential A z is marked byū)[XX] = [ū;Ī c ]. In order to avoid differentiations associated with the Laplace variable s, which may introduce large errors, the entire equation is additionally divided by s N f ilter −1 [2,32]. This corresponds to the division of the unknowns [XX] by different degrees of the variable s. This means that the integrals will now be unknown: From the vector potential udt = Adt and the integral from the current, I c dt, as well as integrals of the higher degrees of these unknowns.
where the number of sub matrices with n 1 rows is equal to the highest degree of the variable s in the IIR filter transmittance, i.e., N f ilter − 1.  Figure 8, is shown in Figure 9.   Figure 10 shows the waveforms of the current I c and choke supply voltage Ew. Figure 11 shows the energies of the supply, sheet metal losses, and resistance losses. The correct balance of these energies indicates the correctness of the calculations. The energies accumulated in the leakage inductances in the air and in the winding area were small enough to be negligible here. In order to calculate the transient waveforms of the Ic current and magnetic flux density B, the losses in the laminated sheets are presented with the use of equivalent filters in the Laplace variable s domain. When multiplying the system of equations of the finite element method by the IIR filter denominator representing the magnetic reluctance of sheets in the domain of variable s, the system of equations became global [39] and required extension of the state variables by their integrals (20) during integration with the Runge-Kutta method [40]. On the other hand, after performing the main FEM calculations for the post-process calculation of the magnetic field strength H in the finite elements, the IIR filter was used in the domain of the variable from the transformation Z. Calculations of this magnitude were, therefore, made locally for each finite element separately. The balancing of these quantities, as illustrated by the zero energy balance, is the evidence of the used method correctness. Similarly to energies, the momentary powers presented in Figure 12 were also balanced. Figure 11. The waveform of the power supply, energy loss in iron, energy loss of resistance, and a zero energy balance for angular frequencies ω = 3141.6 rad s . Figure 12. The momentary power waveform in leakage inductance (L sr ), power in the air area (air), power in the winding area (winding), which all increased 100 times, and the power loss in laminated sheets for angular frequencies ω = 3141.6 rad s .

Conclusions
The issue of appropriate modeling of the phenomena connected with the generation of losses in laminated magnetic cores of electric machines is significant from the perspective of the simulation result correctness of their transient states. Converter systems are being used increasingly more often to supply these machines. Non-sinusoidal character of the supply voltage waveforms additionally increase the power of losses, making the consideration of their influence on waveforms of machine dynamic states even more important, and making the modeling more difficult. A majority of methods for modelling losses in laminated cores that are known in the reference works, are of estimative character and require a previous determination of the waveforms of the machine dynamic states. Due to this fact, they do not provide the possibility to correctly consider the influence of the presence of losses in the core on the machine dynamic properties.
Few methods that have been proposed up to now provide the possibility to consider losses in laminated cores directly during the calculations of machine dynamics and they are characterized by very high computational complexity, which significantly limits their application in practice. This elaboration presents a model of equivalent magnetic permeability of a laminated magnetic core, in the form of an IIR filter. The model was then incorporated in the procedure of fast computations of electric machine dynamic waveforms, by means of the Finite Element Method. With the example of a choke, a series of simulation tests were performed. Correctness verification of the obtained results was done by means of a detailed balance of energy and power in the investigated system.
In compliance with the expectations, the results indicate a strong influence of the skin phenomenon on the value power of losses in the core. This effect is noticeable, in particular, by significant values of the magnetization frequency and a high magnetic permeability of the core. Additionally, these observations prove that the developed model allows for correct reflection of the phenomena that take place in the laminated magnetic core in the calculations of machine dynamics. The acquired values of the power of losses are adequate to the ones included in norms DIN EN 10106 and DIN EN 10207, which describe the magnetic properties of transformer sheets and the results of measurement experiments concerning determination of the power of losses in laminated cores.
A drawback of the presented method is the problematic necessity of global influence of the IIR filter on the system of equations of the Finite Element Method. In the follow-up studies, a solution of this problem is anticipated, by means of developing a model of equivalent magnetic permeability of a laminated core, in the form of a digital filter, possibly to implement in every finite element of the FE model separately. The aim of these efforts is to increase the efficiency of the numerical simulation procedure, so that it might find practical application in analyses of the dynamics of electric machine models.