Iron-Loss and Magnetization Dynamics in Non-Oriented Electrical Steel: 1-D Excitations Up to High Frequencies

This paper presents a thorough introduction to and application of the Parametric Magneto-Dynamic (PMD) model of soft magnetic steel sheets. The theoretical background is reviewed, and two different ways are discussed to account for the viscosity-like effects originating from microscopic eddy currents. This is followed by the theoretical calculation of magnetization dynamics and dynamic hysteresis loops in Non-Oriented (NO) electrical steel. Both classical and viscosity-extended approaches are discussed, with respect to the ability of replicating the dynamic hysteresis loop shape and iron-loss under sinusoidal excitation waveforms up to high excitation frequencies. Comparisons against measurements are analyzed for M400-50A and M235-35A NO electrical steel over a wide range of magnetic flux density and excitation frequencies. The proven accuracy and efficiency of the PMD model make it a valuable tool for the calculation of iron losses in electrical machines and transformers, as well as for an in-depth study of magnetization dynamics in individual laminations.


I. INTRODUCTION
Accurate determination of the iron-loss in the magnetic cores of electrical machines and transformers is of great importance, specifically to allow one to reduce the iron-loss by appropriate design measures in the magnetic circuit, and by the choice of suitable materials [1]- [8]. However, the multifarious and intricate physical mechanisms in different magnetic materials, combined with complex flux waveforms in electrical machines, hamper the development of appropriate methods for predicting the iron-loss and the magnetic field distribution [6], [9]- [12]. This is exacerbated by the fact that material data, i.e., iron-loss and magnetization curves, supplied by the data sheet or acquired through standardized measurements under sinusoidal magnetic flux density waveform, do not apply for the magnetic components of inverter-driven electrical energy converters and rotating electrical machines.
The associate editor coordinating the review of this manuscript and approving it for publication was Bing Li .
Useful engineering methodologies for iron-loss predictions in magnetic laminations, where broad frequency ranges are taken into account, either require a large amount of measured data for parameter identification or are prone to large prediction errors [9]- [26]. The more accurate models are, in general, limited to a specific material (e.g. Non-Oriented (NO) or Grain Oriented (GO) or amorphous electrical steels [27]- [29]), the geometry of magnetic core (e.g. tape-wound or stacked lamination or powder core [30]- [37]) and spatial excitation conditions (pulsating or rotational magnetization). A large quantity of different models and methods, combined with poor knowledge of their limitations, often lead to misuse of these models and inadequate results.
The aim of this work is to present a thorough introduction to the physical background, adequate modeling approaches and their application when analyzing power-loss and magnetization dynamics in NO electrical steels that are subjected to one-dimensional (1-D), i.e., pulsating magnetization. In order to predict the iron-loss and magnetization behavior under such conditions accurately, it is indispensable to utilize a hysteresis model that is coupled directly to an eddy-current (lamination) model. This can be solved by using various approaches, where the lamination model is, e.g., based on Maxwell equations [9], [11], [23], [24], [27]- [36], [38]- [57], saturation wavefronts [58]- [60], equivalent circuits [37], [48], [61]- [63] or is included in the hysteresis model implicitly [43], [64], [65]. A comparison of different 1-D lamination models in terms of mathematical structure, implementation, computational performance, accuracy and spatial discretization can be found in [48]. In general, all the models converge to the same results, depending on the spatial and temporal discretization. However, the computational performance and the flexibility for coupling to external circuits, inclusion into Finite-Element (FE) simulations of electrical machines or simulating voltage-and current-driven regimes varies between the approaches.
It is important to note that all pure 1-D lamination models can reasonably be applied only for devices with predominantly pulsating, i.e., unidirectional, magnetic fields. In many electromagnetic devices (e.g. rotating electrical machines as well as transformers), both unidirectional and rotating magnetic fields are generated. In such cases the direct application of pure 1-D lamination models, i.e., models that take into account alternating magnetization, is inadequate. The correct estimation of rotational losses in all core regions based purely on alternating loss data tends to be impossible [25]. However, these 1-D models are of great importance, because they represent the basis for the synthesis of extended 2-D lamination models that can describe rotational magnetization adequately. A 2-D lamination model can, in general, be obtained by coupling two 1-D models utilizing an adequate vector hysteresis model [9], [30], [31], [44], [66]. Contemporary electromagnetic modeling of complex devices is commonly done by using FE techniques. The discussed extended lamination models can be used either traditionally in the post-processing stage, or can be incorporated into the FE model directly [9], [12], [26], [30], [34], [41], [45], [67]- [70]. A direct incorporation of magneto-dynamic models into the field solution during FE processing allows one to ensure a high accuracy of power-loss and performance calculations of the devices, where the coupling between FE and 2-D lamination models is commonly handled by means of nested iterative schemes. The resulting coupled model, however, increases the computation time and is vulnerable to convergence problems, since it involves two iterative solutions of two coupled nonlinear problems [9], [12], [30], [34], [44], [67]. Therefore, the incorporation of the discussed highly non-linear lamination models, which include vector hysteresis, classical eddy-current and excess effects, into the FE models of electromagnetic devices remains a challenge, and has not yet matured to the level of a standard computational technique. This topic is out of the scope of this paper. As a trade-off between accuracy, robustness, and speed, the power-loss estimation in magnetic cores is, therefore, often applied during the post-processing stage [9], [70], [71].
In this paper, a strongly coupled magneto-dynamic lamination model is presented, that is based on the well-known 1-D Maxwell diffusion equation. The model addresses eddy-current and hysteresis phenomena simultaneously, in such a way that the interplay between skin effect, i.e., macroscopic eddy currents across laminations and hysteresis, can be resolved accurately. This is the so-called Parametric Magneto-Dynamic (PMD) model. The in-depth presentation of this modeling approach is completely new, and offers a different and unique analytical insight into the dynamic magnetization and power-loss mechanisms inside NO steel laminations. The presented model offers the flexibility to implement various inverse hysteresis models, and to be extended to account for the effect of local (microscopic) eddy currents [40]. Further on, the PMD model allows one to simulate current-, voltage-and field-driven problems efficiently, and is ready for coupling to external circuits as well as for use in FE simulations of adequate electrical energy transducers. The presented 1-D modeling approach also represents the basis for modeling of rotating magnetic fields, which is, however, out of the scope of this paper.
In the framework of presented research work, two PMD model versions are introduced, with different levels of complexity. The first is applicable for NO materials with a very fine grain size (classical PMD model). The second accounts, additionally, for the effect of local eddy currents (so-called excess losses) in NO materials with coarse grain size by means of magnetic viscosity (viscosity-extended PMD model). Both PMD model versions are validated, where their predictive power is analyzed by means of comparisons to measured data of two industrial NO electrical steel grades, classified as M235-35A and M400-50A. Further on, an in-depth analysis is presented of the spatialand temporal-distribution of electromagnetic field variables in the electrical steel sheet. The presented theoretical background and analysis allows one to understand the magnetization dynamics, and derive possible measures to reduce the application-dependent iron-loss in NO electrical steels. A comparison of the obtained results for the two materials provides insight into the effect of material parameters (thickness, electrical conductivity, grain structure, . . . ) on the individual power-loss components for 1-D, i.e., pulsating magnetic fields. This paper is organized as follows. In the first section, Section I, a brief introduction is given, where the motivation, focus, goals and structure of the presented research work are discussed. In the second section, Section II, the theoretical background regarding NO electrical steels is presented in great detail, where the physical properties of such materials, fundamental physical background, adequate simplifications and proper boundary conditions are discussed. In the third section, Section III, the PMD model is introduced and discussed thoroughly, where the relations between spatially VOLUME 8, 2020 distributed electromagnetic variables are derived and discussed by means of a spatial discretization of the 1-D Maxwell diffusion equation. In this section, the implementation of non-linear or hysteretic magnetic properties, accounting for additional microscopic phenomena (excess loss) is discussed, considering heterogeneous material properties, power-loss calculation and coupling of the PMD model with external electric circuits. The obtained results are presented in the fourth section, Section IV. Two distinctly different sample sets of NO electrical steel grades and the used measurement setup are presented first. Then, the implementation and identification are presented of both discussed PMD model versions (with and without local, microscopic eddy currents). Next, the iron-loss prediction of the viscosity-extended PMD model is evaluated, where individual power-loss components (hysteresis, Foucault eddy currents and excess eddy currents) are analyzed. Finally, the dynamic magnetization, i.e., the spatial-and temporal-distribution of electromagnetic field variables inside steel sheets are analyzed in the case of non-linear skin effect for the unsaturated, as well as, saturated cases. In the fifth, final section, Section V, a conclusion is given regarding the presented research work.

II. THEORETICAL BACKGROUND
NO electrical steel is used commonly as conductor for magnetic flux in many contemporary electromagnetic devices. The main reason for the widespread use are its favorable price and macroscopic magnetic property, namely magnetic permeability µ, which enables the generation of adequate magnetic fields in such devices [8], [72]- [75]. The magnetic permeability µ of the material is, in general, i.e., for anisotropic magnetic materials, a second rank tensor that represents the magnetic constitutive relation of the material where H represents the vector of the magnetic field and B the vector of the magnetic flux density. The constitutive relation (1) is highly non-linear and hysteretic, i.e., µ(H) and, consequently B(H), depend heavily on the state of the magnetic field H inside the material and its history [72], [74], [76]. The shape of magnetic hysteresis loop reflects a macroscopic, integral picture of complex microscopic magnetization processes inside the soft magnetic material [74], [76]. However, the spatially random domain structure and the small grain size of NO electrical steel allow one to treat µ as homogeneous and isotropic [27], [77]. In the case of a 1-D excitation, the magnetic permeability can, consequently, be described by the scalar function µ. The magnetic permeability µ is favorably high until the material saturates, hence, magnetic materials in devices generally operate in states below material saturation.
The drawback of the iron-based NO electrical steel is the good electric conductivity σ of iron, which, in general, is closer to the conductivity of conductor materials (e.g. copper) than to the conductivity of insulation materials [72], [74]. It reflects complex microscopic conduction mechanics on a macroscopic scale and determines the second constitutive relation of the material where j is the vector of the electric current density and E is the vector of the electric field inside the observed material. In contrast to the magnetic permeability, the electric conductivity σ can often be considered independent of the electromagnetic variables, but is also homogeneous and isotropic. Magnetic cores of many devices are excited with magnetic fields that are changing with time. The change of these variables is, in the discussed electrical steels, adequately low, therefore the quasi-static approximation of Maxwell's equations and represent appropriately the relationships between electric and magnetic subsystems inside the discussed materials [72]- [74]. The presented laws of Faraday (3) and Ampere (4), along with the constitutive relationships (1) and (2), fully describe the electromagnetic phenomena inside the discussed materials. On a macroscopic scale, the electromagnetic phenomena inside NO electrical steels can be basically described as follows: Magnetic cores are excited with electric current in excitation windings. These generate a magnetic field H and, according to the constitutive relation (1), a magnetic flux density B inside the magnetic material. If H and B are changing with time, an electric field E is induced according to Faraday's Law (3). Furthermore, E generates an eddy-current density j due to the good electric conductivity σ of the material (2). These eddy currents generate, according to Ampere's Law (4), a feedback magnetic field, which counteracts the magnetic field inside the material. In this way, both electric and magnetic subsystems are strongly interdependent, where circular eddy currents around the magnetic field represent the unwanted effects that generate additional power loss during the magnetization processes. At higher magnetization dynamics a non-linear skin effect of the magnetic field is also observed inside the material [23], [27], [35], [39], [44]- [46], [59], [78], [79].

A. ONE-DIMENSIONAL MODELING APPROACH
The parasitic eddy currents are generally reduced by two measures. The first is the reduction of the electrical conductivity σ by adding silicon to the iron-based materials [73], [80]. This technique can reduce σ substantially, but, on the other hand, deteriorates the magnetic and mechanical properties of the material and increases the production costs [81], [82]. Therefore, silicon is added only up to a reasonable amount, i.e., for NO electrical steel, a maximum of 3.5 mass-percent is added typically. The second measure is to limit the induced FIGURE 1. A thin, wide and long lamination with its reference Cartesian coordinate system in the center of the cross section of the sheet. eddy currents by reducing the thickness of individual laminations. Usually values between 0.1 mm and 0.5 mm can be found in applications. This comes with the drawback of increased production cost and material-property deterioration [73], [80], [83]. Both approaches are generally used simultaneously in all contemporary magnetic cores that are subjected to time-varying magnetic fields. Focusing on such laminated NO materials, the discussed description of the hard coupled phenomena (1) to (4) can be simplified significantly, as presented below. Figure 1 shows a thin, wide and long NO electrical steel sheet, where the Cartesian reference coordinate system is conveniently defined in the geometric center of the sheet due to the circular nature of the electromagnetic phenomena. The length l m and the width w m of the sheet are far greater than its thickness b (l m b; w m b), where the discussed electromagnetic problem is defined by (1) to (4) only inside the volume of the sheet, i.e., x ∈ − b 2 , b 2 , y ∈ − w m 2 , w m 2 and z ∈ − l m 2 , l m 2 . The lamination of magnetic cores is, ideally, carried out along the direction of the magnetic field, so that the magnetic field is generated perpendicular to the direction normal to the sheet surface [31], [35], i.e., perpendicular to the x-axis in Fig. 1.
The specific geometry of the sheets, combined with the homogeneous and isotropic properties of NO materials, allows one to simplify the diffusion problem substantially [35], [43], [74], [84]. In the case of 1-D excitation, we can assume that the magnetic flux density B inside such a sheet is oriented only in one spatial direction that is perpendicular to the x-axis, e.g., in the direction parallel to the z-axis [31]. Considering that magnetic cores generally represent closed magnetic circuits and neglecting the leakage magnetic field, B can be regarded as independent of the position z, and is fully represented by B = [0, 0, B z (x, y, t)] T . It is worth noting that the observed vector field depends on the positions x and y, as the electromagnetic field is distributed across the x-y cross-section of the sheet. In addition, all the variables can be arbitrary functions of time t. Based on the presented assumptions, Faraday's Law (3) simplifies to The induced electric field E is generated tangentially around the magnetized z-axis, and has, as a result, only components E x (x, y, t) and E y (x, y, t). These are directed along the x-and y-axes, respectively. According to the conductive constitutive relation (2), these two components generate an eddy-current density j, which is flowing circularly with respect to the z-axis. Because w m b, the induced electric-field component E x (x, y, t) [and, consequently, the eddy-current density j x (x, y, t)] is negligibly small compared to the component E y (x, y, t), and can be neglected [35], [36], [45], [74]. Consequently, the second term on the left hand side of (5) can be omitted, as indicated by the strikethrough. This also depicts the main benefit of using laminations inside magnetic cores, as, in this way, the eddy-current power loss is significantly reduced [73]. Due to this simplification, the distribution of electromagnetic variables inside the sheet becomes independent of the y-axis. In other words, it is assumed that induced eddy currents flow only parallel to the y-z-plane, where j = 0, j y (x, t) , 0 T .
According to Ampere's Law (4), the induced eddy currents j y (x, t) generate a circular magnetic field¨¨¨∂ where H x (x, z, t) and H z (x, z, t) are the components of the magnetic field in the x-and z-directions, respectively. Because of the far greater length of the sheet in the z-direction compared to its thickness (l m b), the component H x originating from eddy currents is generated only close to the starting and ending edges of the observed sheet. Therefore, this component can be neglected (or, in the case of a closed magnetic circuit, it is canceled out). Consequently, the first term on the left hand side of (6) can be omitted, as indicated by the strikethrough. In this way, the distribution of electromagnetic variables inside steel sheets becomes independent of the z-axis. This simplification is analogous to neglecting the induced E x on both sides of the strip. Both discussed simplifications are often referred to as neglecting the edge effects [35], [84], whereas the observed sheet is considered infinitely long and wide with respect to its thickness.
Combining the simplified Maxwell equations (5) and (6) by using the conductive constitutive relation (2), the so-called 1-D diffusion equation is obtained [27], [40], [60], [74]. It is important to note that (7) is written in the so-called mixed form, where the functions B z (x, t) and H z (x, t) are connected via the magnetic constitutive relation of the material (1). Both of those functions can be obtained by solving the presented 1-D diffusion equation, where their spatial distributions in one time instant are depicted schematically in Fig. 2 b). It is important to note that functions B z (x, t) and H z (x, t), in general, exhibit different shapes due to non-linearity and hysteresis. From the obtained solution, furthermore functions E y (x, t) and j y (x, t) can be calculated by using (5) or (6). In contrast to B z (x, t) and H z (x, t), the relation between E y (x, t) and j y (x, t) is linear. Therefore, these functions always exhibit the same shape, as depicted in Fig. 2 a). The vector representation of all discussed electromagnetic variables is shown in Fig. 2 c). The term 1-D non-linear diffusion equation implies that all the electromagnetic variables (not only H z and B z , but also E y and j y ) are generated individually only in one (but not necessarily the same) spatial direction. Consequently, a nonlinear diffusion is taken into account, also only in one spatial direction, where the distribution of all variables changes only with respect to the x-axis, as shown in Fig. 2.

B. BOUNDARY AND INITIAL CONDITIONS
The obtained 1-D diffusion equation (7) is a non-linear Partial Differential Equation (PDE). Its solution depends on the geometrical variable x and time t. It describes the electromagnetic field distribution inside a soft magnetic lamination with respect to the position x at different time instants. Adequate Boundary Conditions (BC) at both surfaces of the sheet x = − b 2 and x = b 2 , as well as the Initial Condition (IC) of the magnetic field H z (x, t = 0) have to be determined in order to obtain a proper solution of the problem [39], [42], [43].

1) DIRICHLET BOUNDARY CONDITIONS
Proper BCs are obtained by analyzing the excitation of magnetic cores inside electromagnetic devices. The excitation is realized by energizing one or more windings, which are wound around the magnetic core. According to Ampere's Law, this generates a magnetic field H sur (t), which magnetizes the laminated magnetic core. If a core is assembled from thin sheets that are insulated from each other electrically, the non-linear diffusion phenomena described by (7) are affecting the magnetic field only inside individual laminations [33]. Therefore, the magnetic field on all the surfaces of individual sheets inside the core corresponds to H sur (t) generated by the excitation windings. As a result of the surface magnetic field H sur (t), an electromagnetic wave starts to penetrate the sheet according to (7). In this way, the so-called Dirichlet (or first-type) BCs at both relevant surfaces parallel to the y-z-plane of all the sheets are defined by as depicted in Figs. 2 b) and c) [39], [42], [43]. The symmetry of the Dirichlet BCs at the lower and upper surfaces, combined with the homogeneous electromagnetic properties inside individual sheets, results in a symmetrical distribution of the electromagnetic field inside these sheets [42].
That is why H z (x, t) and B z (x, t) are functions with spatial symmetry with respect to the y-z-plane at x = 0, as shown in Fig. 2 b) and c). In contrast to this, E y (x, t) and j y (x, t) are spatially symmetrical about the sheet's longitudinal z-axis, as presented in Figs. 2 a) and c). This is because E y (x, t) and j y (x, t) are oriented in opposite directions in the upper and lower halves of the sheet, due to the circular nature of the electromagnetic field, as depicted in Figs. 2 c).

2) NEUMANN BOUNDARY CONDITIONS
By integrating (5) over the sheet's thickness b, it can be shown that the difference of the electric field E y (t) over the whole sheet [ Fig. 2 a)] is related to the rate of change of the average magnetic flux density B a (t). This difference is The obtained E y (t) corresponds at the same time to the feedback influence of the changing electromagnetic field, and relates to the induced voltage in all the excitation windings [42], [43].
The discussed E y (t), combined with the odd spatial symmetry of E y (x, t) and j y (x, t), is exploited to obtain proper Neumann (or second-type) BCs [42], [43]. Boundary values due to such symmetry are and where, at the center of the sheets, no electric field is induced, i.e., E y (x = 0, t) = 0. Finally, Neumanns' BCs are obtained 4572 VOLUME 8, 2020 by combining (10) and (11) with (6) and (2). The BC at the lower surface of the sheet is defined by whereas the BC at the upper surface is defined by The discussed spatial symmetry is, furthermore, exploited, to restrict the area of interest to only one half of the sheet thickness (e.g. x ∈ 0, b 2 ) [42], [43]. By doing this, a proper Neumann BC at position x = 0 is determined by 3

) INITIAL CONDITION
Along with the defined BCs, adequate IC is required to solve (7). This can be defined by some function H z (x, t = 0). To simplify the presented analysis, trivial IC can be assumed, where H z (x, t = 0) = 0. Such IC corresponds to the ideally demagnetized state of the soft magnetic material [42], [43].

III. PARAMETRIC MAGNETO-DYNAMIC MODEL
The presented PDE (7), combined with adequate BCs, can be solved for known ICs by using different approaches, such as by applying finite differences (FDs) [28], [42], [43], [49], [51], finite elements (FEs) [9], [34], [42]- [44], [57], mesh-free techniques [30], [36], [41], [48], spatial averaging methods [48], [50], [84], equivalent circuit approximation [62], fractional derivatives [32], [54], or by tackling the problem analytically [46]. Most of these approaches rely on some kind of discretization, whereas, usually, a system of non-linear Ordinary Differential Equations (ODEs) is obtained from the PDE (7). The size of the obtained ODE system depends on the applied method, as well as on excitation dynamics [48]. In this paper, a method is presented where the distribution of magnetic variables inside the sheet is approximated by using piecewise constant functions. The so-called Parametric Magneto-Dynamic (PMD) model of NO steel sheets is obtained using this method [48], [50], [84]. The presented method is based on (5) and (6), and gives a clear understanding of the underlying physical phenomena of non-linear diffusion of the electromagnetic field inside the discussed NO steel sheets. Because all the electromagnetic variables depend on time t, this dependence is, for the sake of clarity in the following equations, not displayed explicitly anymore. Furthermore, as discussed in Subsection II-A, the individual variables are oriented only in one spatial direction. Consequently, all the subscripts that denote spatial directions of individual vari-

A. INFLUENCE OF THE MAGNETIC FIELD ON THE ELECTRIC FIELD
Lets assume a function B (x) that describes the distribution of magnetic flux density across a sheet's half thickness at one time instant, as shown in Fig. 3 a). The distribution B (x) is, in general, changing with respect to time t, which is described by corresponding function ∂B(x) ∂t . This function is shown in Fig. 3 b). According to (5), the distribution E (x) inside the sheet can, in general, be obtained by integrating the function ∂B(x) ∂t from the center of the sheet (x = 0) to an arbitrary position x by where x ∈ 0, b 2 and ξ is a dummy integration variable. In (15) the first BC E (x = 0) = 0 is taken into account, as described in Subsection II-B.
The functions B (x) and ∂B(x) ∂t are, however, not known a priori, because the distribution B (x) is influenced back VOLUME 8, 2020 by E (x). This interdependence can be avoided provisionally if the observed upper half of the sheet is discretized virtually into an adequate number of slices. If the observed sheet is divided equally into N s slices, the thickness b s of each individual slice s is defined by where s ∈ N and 1 ≤ s ≤ N s . If we assume that the first slice (s = 1) starts in the center, and that the last slice (s = N s ) ends at the surface of the sheet, the coordinates of the lower and upper slice borders x ls and x us are defined by and respectively. If these slices are adequately thin, the magnetic flux density inside each slice s can be approximated by its average valueB s , defined bȳ Using (19), B (x) can be approximated by a piecewiseconstant functionB s (x) across the sheet's half thickness. An exemplary case of such a function with N s = 5 slices is presented in Fig. 3 a). Due to the linearity of the derivative, ∂t is approximated by a corresponding piecewise-constant function dB s (x) dt , as presented in Fig. 3 b). Furthermore, it is easy to show that both discussed approximation functions have identical average values B a and dB a dt over the sheet thickness as B(x) and ∂B(x) ∂t , respectively. If the sheet is divided equally in N s slices, e.g., dB a dt inside the sheet can, therefore, be simply calculated by Using the discussed discretization with average values inside individual slices that is defined by (19), allows one to express (15) by where the sum of integrals over all the inner slices [second term on the right-hand side of (21)], as well as an adequate part in the observed slice s [third term on the right-hand side of (21)], is added to the BC in the center of the sheet [first term on the right-hand side of (21)].
The main benefit of the proposed approach is that the spatial integration in (21) becomes independent of the function B (x) and depends only on the spatial discretization of the sheet. This is possible because the values are assumed to be constant across all the slices. By evaluating the integrals in (21), the function E s (x) can be expressed by where a detailed derivation is presented in Appendix-A. The obtained function E s (x) in (22) approximates the spatial distribution of the induced electric field, where this field, inside individual slices E s (x) [second term on the right-hand side of (22)], changes linearly with respect to the position x.
The corresponding distribution E s (x) for depicted dB s (x) dt in Fig. 3 b) is shown in Fig. 3 c). Furthermore, the total change of the induced electric field inside an individual slice E s depends only on the b s of the observed slice and dB s dt inside this slice. This is analogous to the change of the electric field across the whole sheet described by (9).
Finally, the obtained distribution E s (x) satisfies the required BC at both boundaries of the model, as discussed in Subsection II-B. The value of E s (x) in the center of the sheet is always E s (x = 0) = 0, whereas the value at the sheet surface is obtained by evaluating (22) which complies with BC (11). The distribution of the eddy-current density j s (x) inside an arbitrary slice s is obtained by application of the conductive constitutive relation (2) to (22), and is described by The obtained distribution j s (x) is due to the linearity of the conductive constitutive relation (2) also a piecewise-linear function, as shown in Fig. 4

a).
Based on (23), in combination with the constitutive relation (2) and Ampere's Law (6), the change of the magnetic field with respect to the position x at the surface of the sheet is determined by This is shown in Fig. 4 b). A comparison of (25) and (13) shows that (24) always satisfies the prescribed Neumann BCs in the observed sheet.

B. FEEDBACK INFLUENCE OF THE INDUCED ELECTRIC FIELD ON THE MAGNETIC FIELD
The feedback influence of the induced electric field on the magnetic field can be obtained by the application of Ampere's Law (6) to the eddy-current distribution (24), where proper Dirichlet's BCs (8) are taken into account. Because the presented modeling approach is based on taking into account half of the sheet, only the Dirichlet BC at the surface is defined. 4574 VOLUME 8, 2020 This BC is equal to the magnetic field H sur that is imposed by the excitation windings, as discussed in Subsection II-B.
In contrast to this, the value of H (x = 0) at the center of the sheet is unknown. Therefore, the approximated distribution of the magnetic field H s (x) across the sheet can be determined starting from the sheet's surface by where the sum of integrals over all the outer slices [second term on the right-hand side of (26)], as well as an adequate part in the observed slice s [third term on the right-hand side of (26)], is added to the BC at the surface of the sheet [first term on the right-hand side of (26)]. By applying (24) and evaluating the integrals in (26), the function H s (x) can be expressed by A detailed derivation of (27) is presented in Appendix-B. The obtained distribution H s (x) is for the exemplary case with N s = 5 slices presented in Fig. 4 b).
By eliminating j s (x) from (26) and obtaining (27), it is visible that H s (x) depends on dB s dt in all the slices. Furthermore, dB s dt in all the inner slices excites the magnetic field at position x in the slice multiple times [second term on the right-hand side of (27)], whereas the impact dB s dt in the outer slices decreases with respect to the decreasing distance to the sheet's surface [third term on the right-hand side of (27)]. Finally, as presented in Fig. 4 b), the approximated distribution (27) complies with the Dirichlet BC, as well as with all Neumann BCs that are presented in Subsection II-B, which can be proven easily by differentiating (27) with respect to x.

C. DISCRETIZED DIFFUSION EQUATION
The interdependence between the electric and magnetic fields can be described completely by linking the obtained distribution H s (x) described by (27) to the initially assumed dis-tributionB s (x) by using the adequate magnetic constitutive relationship (1). The two discussed functions to be linked have, however, distinctively different properties: H s (x) is a piecewise-square function [ Fig. 4 b)], whereasB s (x) is described by a piecewise-constant function [ Fig. 3 a)]. In order to link H s (x) adequately to the average valueB s across individual slices, the average values of the magnetic fieldH s inside individual slices have to be determined. Based on (27), the average valuesH s are, in analogy to (19), calculated byH In this way, a piecewise-constant distribution functionH s (x) is obtained, as depicted in Fig. 4 b). Because the first four terms on the right hand side of (27) do not depend on position x, the averaging integration of (28) is not affecting those terms. By evaluating the average values of the last two terms on the right hand side of (27),H s (x) can, finally, be expressed byH A detailed derivation of (29) is presented in Appendix-C. Equation (29) not only contains the influence of the induced electric field on the magnetic field distribution and satisfies all the BCs, but enables a straightforward implementation of a non-linear or hysteretic constitutive relationship (1). Because it is written in the mixed form containing both variablesH s andB s , the interdependence is completed when an adequate relationshipH s B s is applied for each slice s. By applying such a relationship, all the valuesB s of VOLUME 8, 2020 the initially assumed piecewise-constant function are related adequately toH s .
Taking into account (29) for each individual slice s ∈ [1, N s ] results in a system of N s non-linear ODEs in mixed form. This system is called the PMD model. The obtained PMD model represents a discretized approximation of the diffusion equation (7), where proper BCs are considered. The PMD model can be expressed in matrix form bȳ The coupling matrix K is symmetric, where its size depends only on the discretization; an N s × N s matrix is obtained by dividing the sheet into N s slices. Furthermore, it contains the geometric and conductive properties of the used material, i.e., thickness b and specific conductivity σ of the sheet. The presented PMD model enables the calculation of the instantaneous non-linear eddy-current diffusion inside an NO steel sheet for arbitrary excitation waveforms [40]. In addition, the obtained PMD model offers interesting insight into the discussed diffusion phenomena. For example, the impact of the rates of change across inner slices on the magnetic field distribution is higher [i.e. is multiplied by a factor (N s − s)], whereas the impact of rates of change across the outer slices decreases linearly to the edge of the sheet. The most interesting term in (29) is, however, the last term on the right hand side, which is also reflected in all the diagonal elements of the coupling matrix K. This term corresponds to the classical eddy-current term 1 3 σ b 2 2 dB a dt = 1 12 σ b 2 dB a dt that is valid for adequately thin sheets where the skin effect of the magnetic field is neglected [60], [74], [78], [85]. Furthermore, if the sheet is approximated by only one slice (N s = 1), the PMD model simplifies to the classical thin-sheet model. If, on the other hand, the skin effect cannot be neglected, the sheet should be divided into several slices [29]. Although the skin effect in individual slices of the PMD model is not taken into account, the three coupling terms [second, third and fourth terms in (29)], describe the skin effect between the slices of the sheet. Therefore, the obtained PMD model can be considered as an extension of the classical thin-sheet model, where the basic principles of the thin-sheet model are applied on all the virtual slices of the PMD model, but, furthermore, the skin effect inside the sheet is taken into account.

D. APPLYING NON-LINEAR AND HYSTERETIC PROPERTIES
The non-linear relationshipH s B s in individual slices can be taken into account by applying a proper non-linear function or a hysteresis model [34], [38], [40], [41], [78], [79], [86], [87]. The mixed form of the ODE system (30) offers a straightforward implementation of primal, i.e., inversely formulated hysteresis models, where H is calculated from the known B or dB dt [88]. Over time, lots of different hysteresis models were developed because magnetic hysteresis represents the most complex physical phenomena inside soft magnetic materials [74]- [76]. These range from complex physical-based to simpler, completely phenomenological hysteresis models [40], [77], [86]- [92]. For application in the PMD model, the inverse hysteresis model should replicate the so-called static hysteresis curves properly [93].
The choice of a suitable hysteresis model to be incorporated into the PMD model depends, in most cases, on the application of the coupled lamination model. For example, simpler hysteresis models can be applied in applications with sinusoidal excitations. These are models that cannot replicate complex magnetization curves [e.g. offset minor loops under pulse-width-modulation (PWM) excitation], but are easy to parametrize and computationally efficient. On the other hand, complex magnetization excitations require the application of more detailed hysteresis models, (e.g. history-dependent models), that can achieve high accuracy but require a more elaborated identification, and are computationally more expensive. A detailed presentation of various suitable hysteresis models is given in [40].

E. ACCOUNTING FOR ADDITIONAL MICROSCOPIC PHENOMENA
The homogeneous and isotropic 1-D lamination model presented in Section II allows one to model the non-linear diffusion of the electromagnetic field inside NO steel sheets efficiently. The assumption of a perfectly homogeneous field distribution, i.e., an infinitely fine domain structure, allows one to neglect the effect of microscopic local eddy currents that are induced around the moving magnetic domain walls inside the material [27], [94]. However, these local or, so-called excess eddy currents, have to be taken into account when considering most real NO materials, in particular when modeling NO steel sheets with a coarser-grained structure [27], [95]- [97].
The macroscopic effect of these micro-scale eddy currents is expressed as an additional lag of the magnetic flux density B behind the magnetic field H , whereas the micro-scale eddy currents locally dampen changes of the field inside the lamination. Various approaches and models, which range from purely phenomenological to more physically based, were developed in order to model such effects. Different models are presented in [51]- [56], [65]. In this paper, a physical-based approach was considered, as the discussed macroscopic effect resembles viscous-like friction. On that account, the PMD can be extended by using the so-called magnetic viscosity [27], [51], [52].
The size range of magnetic domains inside NO steel sheets that are used in contemporary electromagnetic devices is considerably smaller than the thickness of the lamination [75]. Based on the adequately small and randomly distributed magnetic domains, one can assume that such domains do not change nearly as significantly as, e.g., in the case of GO materials. The induced eddy currents around such moving domain walls are, therefore, generated only by the local magnetic flux density, and, consequently, influence back the magnetic field in the same (local) spatial region. Based on this, an additional component of the so-called average viscous magnetic field H vs is added in each slice of the PMD model. Such a field componentH vs can, in general, be described bȳ where δ is a directional variable, g(B s ) represents an adequate function that controls the impact of the viscous magnetic field H vs and α is a model parameter. By introducing an adequate complex function g(B s ) and upgrading the parameter α to be dependent onB s , virtually any dynamic loop shape can be reproduced [28], [29]. In this way, the presented PMD model extension can also be used practically for modeling of other magnetic materials (e.g. GO electrical steels) despite that the underlying physical mechanisms are significantly different [29]. However, by applying such functions, the complexity of the obtained model and its parameter identification is increased significantly, and the physical interpretation of its parameters is lost. Therefore, its practical value is, consequently, reduced. For NO materials, a simple, physical-based function g(B s ) can be applied [27], [29], [57]. In this way, the viscosity termH vs inside individual slices can be described efficiently bȳ where R m is a constant model parameter and B sat is the saturation flux density of the material [28], [29], [57]. In the presented physical-based approach it is considered that magnetic domains change their sizes only when the NO material is not saturated. This is taken into account by the middle term on the right hand side of (33), which is enclosed in parentheses. The structure and mobility of domain walls is reflected in parameter R m . Furthermore, the induced viscous magnetic field depends exponentially in respect to the changes of local magnetic flux density, whereas consistency with the statistical loss theory can be accomplished if α = 2 is selected. Application of (33) allows one to appropriately replicate the lag of the magnetic flux densityB s behind the applied field inside individual slices for NO electrical steels. The comparison between taking viscosity effects into account and neglecting them is for the exemplary case presented in Fig. 5. By incorporating (33) into the PMD model (30), the original system of non-linear ODEs becomes a system of non-linear Differential Algebraic Equations (DAEs). Iterative numerical techniques are required to solve such a system [98], [99]. This can be avoided, if (33) is extended pragmatically to an ODE bȳ where τ v is a time constant. The dynamic term introduced in (34) enables the expression of the viscosity extended PMD model as a system of ODEs. Furthermore, it delays the impact of the magnetic viscosity so that it does not change instantaneously, but continuously with a prescribed dynamic. Implementing magnetic viscosity effects according to (33) or (34) in all the slices and rearranging (30), the viscosity VOLUME 8, 2020 extended PMD can be expressed by where the applied magnetic field H sur equals the sum of three components (the so-called field-addition principle): Magnetic field H h due to static hysteresis, magnetic field H e due to induced macroscopic eddy currents, and magnetic field H v = H v1 ,H v2 , . . . ,H vNs T due to micro-scale local eddy currents.

F. CONSIDERING HETEROGENEOUS MATERIAL PROPERTIES
The spatial discretization of the PMD model enables incorporation of heterogeneous properties inside the sheet. The properties of the material can be described by a piecewise-constant function across the sheet thickness. This can be achieved by assuming different electrical conductivities σ , spatial-dependent relations betweenH s andB s , as well as locally-varying viscous properties inside individual slices in the sheet. Furthermore, the discretization of the model can, theoretically, be arbitrary, where the sheet is not necessarily divided into equally thick slices. According to the distribution of discretization, the elements in the coupling matrix K also change their values. Such heterogeneous properties are, however, still symmetric to the sheet's half thickness, as only half of the sheet is modeled.

G. POWER LOSS CALCULATION
The basic version of the PMD model enables calculation of the so-called classical power-loss components, due to static hysteresis and due to non-local (macroscopic) eddy-currents.
The viscosity extended PMD model additionally enables calculation of the power-loss component due to local (microscopic) eddy currents, which is known as the so-called excess power loss. The instantaneous power-loss components that correspond to all these three phenomena inside the sheet can be determined based on the local electromagnetic variables inside individual slices. In this way, besides the total power loss in the sheet, the spatial distribution is also obtained of power losses across the observed sheet. For the sake of generality, all the power and energy components are, in this paper, regarded as specific values per unit mass, i.e., are determined in W/kg or J/kg, respectively. The instantaneous power required to generate the non-local (macroscopic) eddy currents in slice s can be obtained by integration of the product of the distributions j s (x) and E s (x) over the volume V s = w m b s l m of the observed slice. Because only half of the sheet is taken into account, the obtained power for one slice has to be considered twice. Furthermore, this power is divided by the mass m of the sheet to obtain specific value per unit mass p es by The mass of the sheet can be determined by m = 2 N s V s ρ, where ρ represents the volumetric mass density of the sheet. By considering this and inserting (24) into (36), p es can be determined by The obtained p es with (37) is always positive, therefore, it represents the specific power loss due to non-local eddy currents. Furthermore, p es can be considered as an extension of the classical eddy-current loss of the thin-sheet model [50], [59], [74], [95]. The last term on the right hand side of (37) represents the classical eddy-current loss inside the observed slice s, whereas the first two terms represent the increase of the loss due to the non-linear eddy-current diffusion.
The instantaneous power per unit mass p hs that is required to overcome the static hysteresis in slice s considering both halves of the sheet can be calculated by The static magnetization phenomena of magnetic materials include both reversible, as well as irreversible energy conversions. In the reversible processes, the energy is stored in the magnetic field, and is returned back to the power source (or is later dissipated elsewhere in the electromagnetic system) when magnetization occurs in the reverse direction. In contrast to this, the energy in the irreversible processes is converted into heat. Only this part of the energy is contributing to the static hysteresis power loss. Consequently, the instantaneous power per unit mass p hs can be both positive, as well as negative, where the specific power loss can be calculated only for completed magnetization cycles, as discussed below. Analogous to p hs , furthermore, the instantaneous power per unit mass p vs that is required to generate the local, microscopic eddy currents, can be calculated by The viscosity processes are irreversible, therefore p vs is always positive. It represents the specific power loss due to local eddy currents, often regarded as the excess power loss [72], [74]. The distribution of all the discussed instantaneous power components across the sheet's thickness is, for the exemplary case corresponding to Fig. 3, Fig. 4 and Fig. 5, presented  FIGURE 6. Distribution of power-loss components inside a half sheet's thickness that is divided into N s = 5 slices: a) total power loss p ts and non-local (classical) eddy-current power p es b) hysteresis power loss p hs , and c) local eddy-current loss due to magnetic viscosity p vs .
in Fig. 6. In the presented results, power components for both halves of the sheet are taken into account. The theoretical instantaneous specific powers p es , p hs and p vs can be used to calculate the individual components of energy per unit mass W •s as well as average power loss components P •s . These can be calculated by where t represents an arbitrary time period that starts at t 1 and ends at t 2 and corresponds to a magnetization cycle, and • represents a placeholder for index h, e or v. The total instantaneous powers p e , p h and p v , energies W e , W h and W v , and average powers P e , P h and P v for the entire sheet are simply obtained by summing all the components of the slices.

H. COUPLING WITH EXTERNAL ELECTRIC CIRCUITS
External electric circuits that excite a magnetic core assembled from discussed NO steel sheets are usually described by the electric circuit theory that represents a coarse approximation of the field theory [42], [99], [100]. The circuit theory is based on the two Kirchoff Laws, combined with adequate constitutive relations of lumped circuit elements. The basic elements have two poles, and include resistors, inductors, capacitors and various sources, whereas the behavior of each such element is described by only two scalar quantities, which are the current through and the voltage across the element. A systems of ODEs is usually obtained by modeling the external excitation circuits using the circuit theory. Such system of ODEs can be extended with the ODEs describing the PMD model, whereas the obtained coupled system is solved as a whole.
The presented PMD model can be coupled with an external circuit exploiting the BCs discussed in Subsection II-B. In this way the PMD model is considered as a two-pole element that can be applied into an adequate circuit model. The current i through this element is, according to Ampere's Law, connected to H sur by where N is the number of turns of the excitation wind ing [36], [42]. The voltage across such an element equals, according to Faraday's Law, the induced voltage u i inside the excitation winding. The induced voltage u i is generated by the changing magnetic flux inside the whole magnetic core that is assembled from NO lamination, and is described by where A Fe is the cross section of the magnetic core [42].

IV. RESULTS
In this section, the parameter identification and validation of the PMD model are presented for two distinctly different sets of samples of NO electrical steels. The first material is classified as M400-50A. This material contains 2.4 wt% mass of silicon, and has a nominal total thickness (including the insulation layer) of 0.5 mm. The second material is classified as M235-35A. This material contains 3.5 wt% mass of silicon, and has a nominal total thickness (including the insulation layer) of 0.35 mm. The two electrical steel grades were chosen consciously, for their widespread use in applied engineering and significant difference in their physical properties. The first grade is on the thicker side of the production range and has lower silicon content, which resulted in higher classical eddy current losses. Such a grade is used mainly in industrial applications, and allows one to study the effect of a large amount of non-local eddy currents. The second material is significantly thinner and has higher silicon content, which resulted in lower eddy current losses. Such a material is used commonly in automotive and other high-efficiency applications. Furthermore, the micro-structure of both steels is different. The measured geometric and constitutive properties VOLUME 8, 2020  of the evaluated samples are presented in Table 1 and Table 2, respectively.

A. EXPERIMENTAL EVALUATION
Both materials under study were evaluated individually experimentally by using an Epstein frame within a computer-aided setup in accordance with the International Standard IEC 60404-2. For the experimental evaluation a mixed set of Epstein strips were used, i.e., including half of the strips oriented in a rolling direction, as well as half of the strips oriented in a transversal direction. The samples were evaluated experimentally for a cyclic average magnetic flux density B a (t) = B max sin(2πft), that is controlled to be sinusoidal. In this way, dynamic hysteresis loops for maximum average flux densities B max between B max = 0.5 T and B max = 1.6 T were obtained in the frequency range between f = 5 Hz and f = 1000 Hz. Both measuring ranges were limited, due mainly to the abilities of the measurement system, where the measurements out of these ranges did not fulfill the required form-factors, or were unreliable.

B. IMPLEMENTATION OF THE PMD MODEL
For the detailed model analysis, two versions of the PMD model were implemented by using the commercial simulation software Matlab/Simulink: • the first version is based on (30), and represents the so-called basic PMD model, where only the macroscopic (non-local) eddy currents were considered, • the second version is based on (35), and represents the viscosity-extended PMD model, where both the macroscopic (non-local), as well as microscopic (local) eddy currents inside the sheets are taken into account.

1) APPLIED STATIC HYSTERESIS MODEL
As discussed in Subsection III-D, the mixed forms of both PMD model versions require an inversely formulated hysteresis model for implementation. The ODE-based Tellinen hysteresis model is applied in this paper [79]. The advantages of this hysteresis model are its easy implementation, its computational efficiency and its direct usage of measured major static hysteresis loops. Consequently, it is able to replicate major loops of all shapes. In contrast to this, many other hysteresis models cannot reproduce measured static loops correctly, or become very complex, in order to achieve adequate accuracy. The accuracy of the Tellinen model when reproducing minor loops is comparable to other well-known hysteresis models, such as e.g. the Jiles-Atherton model [40]. The deviation of the symmetric minor loops is also very dependent on the magnetic properties of the considered material, where these deviations are bigger when modeling materials with steeper static major hysteresis loops (higher permeability). The main drawback of the Tellinen model is, however, that it is history-independent, i.e., has no memory. Due to this, complex magnetization curves (generated by, e.g., PWM-excitation waveforms) cannot be reproduced correctly [40].

2) DISCRETIZATION OF THE PMD MODEL
The spatial discretization of the PMD model can, in general, be adapted to the excitation dynamics: If there is no skin effect (i.e. penetration depth is bigger than b 2 ), then N s = 1 is adequate. In this way, the influence of macroscopic eddy currents is taken into account by using the well-known classical eddy-current term 1 12 σ b 2 dB dt [59], [85]. However, when the non-linear skin effect is present, N s has to be increased accordingly. In the case of the basic PMD model approximately 1 slice per penetration depth is required. If the PMD model is extended using magnetic viscosity, the level of required discretization can be reduced significantly. As discussed in Subsection III-E, the introduction of magnetic viscosity dampens quick changes inside individual slices, therefore, a coarser discretization is already adequate [29].

3) NUMERICAL PROPERTIES AND SOLVING
The advantage of the obtained ODE-based PMD models is that the solution can be calculated directly by using the backward differentiation formula for direct integration. In this way, an iterative solving is avoided at each time step. When implementing a hysteresis model into both versions of the PMD model, the obtained coupled models classify as stiff ODE systems, which require the use of adequate numerical methods for accurate and efficient solving [98]. The implemented models were solved by using the Simulink's built-in variable-step solver for stiff ODE systems that is named ode23tb. This solver is an implementation of TR-BDF2, an implicit Runge-Kutta formula with a trapezoidal rule step as its first stage, and a backward differentiation formula of the order two as its second stage [101].

4) COUPLING WITH AN EXTERNAL CIRCUIT
Finally, it is important to note that the discussed models can be implemented either as so-called direct-driven models, or as so-called indirect-driven (or feedback) models. In the first case, the input variable for the model is H sur or i, where the field has no feedback influence on those two variables. In the second case, the input variable for the model is, e.g., voltage u in the excitation winding. This voltage generates a current i in the excitation winding and corresponding H sur on the magnetic core. The resulting electromagnetic field, in turn, generates induced voltage in the excitation winding u i that has a feedback influence on the current i and, consequently, H sur . In this paper, the voltage driven models that are coupled with an excitation winding are implemented, where measured values of u are used as the input variables. In this way, the structure and operation of the measurement system that was used for validation of the results is taken into account.

C. BASIC (CLASSICAL) PMD MODEL 1) PARAMETER IDENTIFICATION
The main advantage of the basic (classical) PMD model (30) is that all its parameters are physically-based, and can be obtained directly from measurements. Consequently, the presented data in Table 1 or Table 2 include all but one physical parameter of the PMD model (30). The missing material property is the non-linear and hysteretically changing magnetic permeability of the material, which is introduced implicitly into the PMD model by applying the discussed Tellinen hysteresis model. To obtain a proper hysteretic relationship for the hysteresis model, the electrical steel samples have to be evaluated experimentally under so-called static magnetization conditions. These conditions are achieved by the application of an extremely low rate of change of the magnetic flux density, where all the discussed dynamic effects caused by macroscopic, as well as microscopic, eddy currents are negligible [93]. Such measurements are very sensitive and are very time consuming. In this paper, the major static hysteresis loop of the steel sheets is determined by extrapolating dynamic hysteresis loops that are measured at low excitation frequencies (i.e. f = 5 Hz and f = 10 Hz) at B max = 1.6 T. This extrapolated major loop is used as a look-up table in the Tellinen hysteresis model.

2) EVALUATION VERSUS MEASUREMENTS
The first version of the coupled PMD model was evaluated based on the obtained parameter set. This version takes into account only the macroscopic eddy currents. The comparison between the theoretical and measured dynamic hysteresis loops at B max = 1.5 T for samples of material M400-50A is shown in Fig. 7.
The results presented in Fig. 7 show that the discussed model predicts dynamic hysteresis loops that have a similar shape as the measured loops over the whole analyzed frequency range, where the calculated dynamic loops are narrower compared to the measured ones. The shape of the dynamic loops at different excitation frequencies is strongly connected to the non-linear eddy-current diffusion, or gradual electromagnetic wave penetration (phase shift of electromagnetic field -see Subsection IV-F), across the sheet thickness. The penetration of the wave is caused by the induced eddy currents and influenced by non-linear magnetic properties [49], [58], [59], [78], [85]. The eddy currents, in general, increase the power loss inside the sheet, which is reflected in the increased surfaces of dynamic loops. At low excitation frequencies, where the skin effect or non-linear eddy-current diffusion is negligible (i.e. penetration depth δ of the magnetic field is bigger than half of the sheet thickness), the inflated dynamic loops have almost parallel branches to the static loop. At higher excitation frequencies (i.e. penetration depth δ of the magnetic field is smaller than half of the sheet thickness), the inflation of the dynamic loops is not parallel to the static loop anymore. In such cases, the inflation of dynamic loops due to phase shift of the electromagnetic field inside the sheet increases with ongoing magnetization in one direction. Therefore, the inclination (as well as the effective magnetic permeability) of the dynamic loops decreases with increased frequency. The saturation practically eliminates the phase shift. As a result, the curves in the saturation region become parallel to the static loop again. The penetration depth δ for the analyzed samples was estimated to δ = b 2 at approximately f = 70 Hz. Consequently, both the theoretical, as well as measured dynamic loops at f = 10 Hz in Fig. 7, are almost parallel to the static hysteresis loop. In contrast to this, all the dynamic loops at higher frequencies in Fig. 7 have magnetization curves that are not parallel to the static loop anymore. This deviation is increasing with increased excitation frequency. The presented results in Fig. 7 show that the non-linear skin effect, or phase shift of the magnetic field inside the sheets due to induced macroscopic eddy currents, is the dominating process that inflates the dynamic loops. However, it is also obvious that the predicted dynamic loops are narrower compared to the measured ones. The difference between the theoretical and measured dynamic loops in Fig. 7 is highlighted by the hatched area between the loops at individual frequencies.
The obtained results show that the deviation of the predicted vs. measured hysteresis loops is almost constant when the sheet is not saturated, where the deviation vanishes when the sheet is saturated. Therefore, this deviation can be linked to the level of B inside the slices, as described by (33). The discussed deviation is, furthermore, frequency dependent, where the error increases exponentially from area 1 to area 4 . This deviation can be attributed to microscopic eddy currents that are induced around the moving domain walls, as described in Subsection III-E by (33). To improve the accuracy of the first version of the PMD model, the magnetic viscosity was added to each individual slice.

D. VISCOSITY-EXTENDED PMD MODEL 1) PARAMETER IDENTIFICATION
Based on the presented results in Subsection IV-C, the simplified extension for magnetic viscosity is described by (34), and represents a proper upgrade of the PMD model. Such an extension requires the identification of four additional parameters, namely B sat , α, R m and τ v . In this paper, B sat was calculated based on the physical properties of the analyzed samples, and parameter α was set to α = 2, in order to keep consistency with the statistical theory of excess loss [39], [95], [102]. The time constant τ v was determined adequately small, so that it did not influence the obtained results, but, at the same time enabled that the model was solved efficiently as a system of ODEs. In this way, only the parameter R m was left to identify. The latter controls the parallel viscosity (excess loss) expansion of the dynamic loops [57]. R m in this paper was determined simply by minimizing the deviation between the total power loss of the calculated and measured dynamic hysteresis loops in one evaluation operating point at f = 200 Hz and B max = 1.5 T. All the identified parameters are, for both discussed materials, gathered in Table 3 and Table 4, respectively.

2) EVALUATION VERSUS MEASUREMENTS
The viscosity-upgraded PMD model version was evaluated using the sets of parameters presented in Table 3 and Table 4. The comparison between theoretical and measured dynamic  hysteresis loops at B max = 1.5 T is presented in Fig. 8 for material M400-50A.
The obtained results show a very good agreement with the measured dynamic hysteresis loops over the whole analyzed frequency range. Small deviations between theoretical and measured magnetic hysteresis loops are apparent, especially in the knee region of the hysteresis curves. Based on a sensitivity analysis, these deviations can be attributed to the accuracy of the static hysteresis description. As the static hysteresis was estimated by extrapolation from the dynamic loops at f = 10 Hz and f = 5 Hz, the biggest uncertainty was introduced in the knee region. The performed sensitivity analysis showed that an accurate static hysteresis model is indispensable for prediction of complex magnetization curves that are generated e.g. when PWM-excitation waveforms are applied.

E. POWER-LOSS CALCULATION 1) EVALUATION VERSUS MEASUREMENTS
In this section, the viscosity-extended version of the PMD model is evaluated versus the total measured power loss P t,meas . The evaluation was performed in the frequency range up to f = 1 kHz and a maximum flux density range from B max = 0.5 T to B max = 1.5 T for both considered materials. The obtained results for material M400-50A are presented in Fig. 9, whereas the results for material M235-35A are presented in Fig. 10.
The comparison between theoretical and measured total power loss (P t and P t,meas , respectively) shows a good agreement for both materials in the whole analyzed frequency range [ Fig. 9 and Fig. 10; subfigures a), e) and i), respectively], as well as over the analyzed range of B max [ Fig. 9 and Fig. 10; subfigures b), f) and j), respectively]. The absolute error ε a shown in Figs. 9 e) and f), as well as Figs. 10 e) and f), was determined in each evaluation point as the difference between measured and theoretical total power loss by whereas the relative error ε r depicted in Figs. 9 i) and j) as well as Figs. 10 i) and j) was determined by The obtained results for M400-50A in Fig. 9 show that the relative error ε r in the analyzed ranges was lower than 5 %, whereas ε r was, in most evaluation points, even lower than 1 %. The biggest deviations occurred at low excitation frequencies, where the static hysteresis phenomena dominated the generation of power loss. This deviation can be attributed to the uncertainty of the used static hysteresis description, as discussed in Subsection IV-C. In comparison to M400-50A, results for M235-35A in Fig. 10 show similar trends, but slightly bigger deviations versus measurements. The relative error ε r in this case was, in both analyzed ranges, mostly lower than 5 %. However, the biggest error ε r was again obtained at low f and B max , and was as high as 17 %. For these evaluation points an additional analysis was carried out of the accuracy of the static minor loops predicted by the Tellinen hysteresis model. The results of this analysis showed that the accuracy of calculated static hysteresis loops at low B max was significantly lower for M235-35A compared to M400-50A. This analysis confirmed again that the obtained deviation in the calculated power loss can be attributed to the uncertainty and poor accuracy at low VOLUME 8, 2020 B max of the used static hysteresis description, as discussed in Subsection IV-C.
The performed evaluation has confirmed that the presented modeling approach is adequate for calculation of power loss in a broad f and B max ranges under 1-D excitation of different NO steel sheets. The performed analysis, furthermore, emphasized the importance of accurate description of the static hysteresis. All the dynamic variables are directly dependent on the static hysteresis description. Therefore, an accurate static hysteresis model is necessary, especially when distorted magnetization of NO steel sheets is considered, where symmetric and offset minor loops are generated.

2) ANALYSIS OF POWER-LOSS COMPONENTS
As presented in Subsection III-G, the theoretical total power loss P t is originating from three distinct loss phenomena: 1) the static hysteresis phenomenon generates the so-called hysteresis power loss P h that is presented in Subfigures c) and d), 2) induced macroscopic eddy currents generate eddycurrent power loss P e that is presented in Subfigures g) and h), and 3) microscopic phenomena generate the viscosity (excess) power loss P v that is presented in Subfigures k) and l).
These power-loss components are presented for M400-50A in Fig. 9, and for M235-35A in Fig. 10, respectively. The obtained results show that material M235-35A has slightly lower P h compared to material M400-50A. Lower P h can be attributed to the narrower static hysteresis of M235-35A due to smaller amount of dislocations. The difference between P v is bigger, whereas lower P v reflects lower microscopic eddy currents and, therefore, the smaller and smoother domain structure of M235-35A. This is also reflected in the viscosity parameter R m . P e represent the biggest difference between both materials. Higher silicon content (resulting in lower σ ) and thinner lamination (lower b) of M235-35A resulted in significantly lower P e . The presented results, furthermore, enabled us to analyze the dependence of individual power-loss components in respect to f and B max . Results for both materials were similar and lead to the same conclusions. The dependence of P h in respect to f can be, based on the representation in Figs. 9 c) or Fig. 10 c), visually described as a linear function. However, this dependence was linear only in the f and B max ranges, where the skin effect is negligible (i.e., at low excitation frequencies or at high B max , where the skin effect is limited by saturation). The presence of non-linear skin effect introduced a slight non-linearity in the discussed dependence, where the skin effect increased P h in comparison to the classical linear relationship. This non-linearity increased when B max decreased, because, in these cases, the skin effect was more pronounced. When analyzing the dependence of P h in respect to B max , it is clear that this dependence was highly non-linear. At low excitation frequencies, this dependence reflected the non-linear increase of the static hysteresis loop surfaces with increasing B max , whereas at higher frequencies, this was additionally influenced by the skin effect. For both analyzed materials this dependence could not be approximated adequately with a simple exponential function using one term, e.g., P h ∝ B n max , where n is a real value exponent. The dependence of P e in respect to f is also influenced by the non-linear skin effect. This dependence is, according to the statistical theory of power loss, described by a quadratic function without the linear term, e.g., P e ∝ f 2 [74]. Similar to discussed dependence of P h , P e can be fitted in accordance with the statistical theory of power loss only when f is adequately low and, consequently, there is no skin effect. A quadratic fit can be applied pragmatically in the whole analyzed frequency range when the skin effect is limited by saturation (e.g. B max = 1.5 T or higher). At lower values of B max this dependence can be approximated over the whole analyzed frequency range with an exponential function P e ∝ f n , where the value of the exponent is decreased (i.e., n < 2). For both materials, the value of the exponent decreased to around n ≈ 1.75 at B max = 0.5 T. At this B max skin effect is significant in the majority of the analyzed frequency range. Consequently, it can be concluded that the presence of skin effect decreases P e compared to the statistical theory. Similar to P h , the dependence of P e in respect to B max is, for both materials, hard to adequately approximate with a single term exponential function over the whole range. In the presented cases, a loose approximation was obtained with P e ∝ B 1.8 max . Lastly, also, dependence of P v with respect to f is influenced by the skin effect. When the skin effect is negligible, the dependence can be approximated by an exponential function P v ∝ f 1.5 , which corresponds to the statistical theory of power loss. In the broader frequency range, a better approximation is found if the value of the exponent is n ≈ 1.6. In contrast to P h and P e , this dependence was found independent of saturation. The dependence of P v in respect to B max was, again for both materials, hard to approximate adequately with a single term exponential function over the whole B max range.
The performed analysis shows the influence of individual material parameters on the different power-loss components. The obtained power loss is in accordance with the statistical loss theory only in excitation cases where the non-linear skin effect is negligible [46].

F. DYNAMIC MAGNETIZATION INSIDE NO ELECTRICAL STEELS
The validated viscosity-extended PMD model was, furthermore, used to evaluate the instantaneous spatial distribution of all the electromagnetic variables inside the discussed steel sheets. The results are presented for samples of M400-50A NO electrical steel. For the purpose of showing the discussed time dependence as well as spatial distributions of all the variables, the steel sheet was over-discretized; N s = 10 slices were used to obtain smoother time and spatial dependences of electromagnetic variables. Especially interesting magnetization behavior was obtained at higher excitation frequencies, where the non-linear skin effect was generated. Therefore, results are shown for f = 400 Hz and two distinct magnetic flux density levels. These are B max = 0.5 T and B max = The first presented case is shown in Fig. 11. In this case, the electromagnetic variables inside the sheet are not limited by saturation as B max = 0.5 T. Significant non-linear skin effect is observed due to the adequately high excitation frequency f = 400 Hz. The magnetic field at the start (t/T p = 0, B a = B max ) and end (t/T p = 0.5, B a = −B max ) of each half period was forced to the edges of the sheet, as shown in Figs. 11 a), b), c) and d). At the center of the sheet (slice 1), the maximum density of magnetic flux was only B max1 = 0.35 T, whereas near the surface (slice 10) B max10 = 0.81 T. Furthermore, a phase shift between B s inside individual slices was observed, where magnetization started at the surface of the sheet and gradually continued into the sheet to the center. This corresponds to an electromagnetic wave that penetrated the observed sheets from all the surfaces to their centers. Furthermore, due to the non-linear and hysteretic magnetic constitutive relationship, the spatial distributions of magnetic flux density B (x/b) shown in Fig. 11 b) and magnetic field H (x/b) shown in Fig. 11 d) have significantly different spatial profiles. These magnetization profiles show that the influence of discussed macroscopic, as well as microscopic phenomena, results in a non-uniform magnetization across the sheet at virtually all time instants.
Both macroscopic and microscopic phenomena are dependent on the rates of change of magnetic flux densities dB s dt inside the sheet, therefore, these phenomena are considered as dynamic. The rates of change of magnetic flux densities dB s dt across the sheets are shown in Figs. 11 e) and f). Because of the skin effect, dB s dt are significantly higher in the slices near the surface compared to those in the center. VOLUME 8, 2020 The resulting distribution of the viscous magnetic field H v is, according to (34), a direct consequence of the change rates dB s dt inside individual slices. As discussed in Section III-E, in the case of NO steel sheets, the size range of magnetic domains is adequately small, therefore, local rates of changes dB s dt generate local microscopic eddy currents that are independent of rates of change in other slices. These microscopic eddy currents generate a local viscous 4586 VOLUME 8, 2020 magnetic fieldH vs that influences back the magnetization across the steel sheet. Consequently the calculated time profiles shown in Fig. 11 g), as well as spatial profiles shown in Fig. 11 h), resemble profiles of dB s dt in Figs. 11 e) and f), respectively. Furthermore, according to (33) or (34), the microscopic phenomena are damping the changes of electromagnetic variables inside the steel sheets. Therefore, a coarser discretization of sheets can also be applied if microscopic phenomena (viscosity) are taken into account in comparison to the case when these effects are neglected [29], as discussed in Subsection IV-B.
In contrast to this, the macroscopic induced electric field E and, consequently, eddy current density j, cannot be calculated based on local rates of changes dB s dt individually. According to (22), E s is generated by dB s dt in the observed slice, as well as the rates of change inside all the inner slices. Furthermore, the induced electric field is oriented in different directions in the upper and lower halves of the sheet, as presented in Figs. 11 i) and j). The induced electric field is, due to the symmetry of the discussed problem, zero in the center of the sheet, and changes non-linearly to both surfaces of the sheet. The presented results show that the use of the classical eddy-current factor 1 12 b 2 dB a dt in the case of non-linear skin effect is not adequate, as, in that case, E(x) is always approximated with a linear function across the steel sheet, as discussed in Subsection III-C.
The eddy-current density j generated by E s is a linearly scaled version of profiles shown Figs. 11 i) and j), therefore, it is not presented explicitly. The induced j generates a magnetic field H e that is presented in Figs. 11 k) and l). Based on (26), H es inside individual slices is generated by j s in the observed slice, as well as by eddy-current densities in all the outer slices. Consequently, in contrast to all other electromagnetic variables, H e is the highest in the center of the sheet, and lowest at both surfaces. By eliminating j s in (27), it was shown that H es depends on dB s dt in all the slices of the sheet. Consequently, the profiles of H e do not resemble individual profiles of dB s dt as in the case of H v , but are significantly different.
Based on (37)-(39), the individual power-loss components were calculated that are required for magnetization of the steel sheet. The instantaneous total input power p t that is presented in Figs. 11 s) and t) consists of three components that correspond to all three distinct phenomena inside the observed sheet: 1) the instantaneous power p e is needed to generate the induced macroscopic eddy currents and is shown in Figs. 11 m) and n), 2) the instantaneous power p h is needed to overcome static hysteretic properties and is shown in Figs. 11 o) and p), and 3) the instantaneous power p v is needed to generate the induced microscopic eddy currents and is shown in Figs. 11 q) and r).
According to (36), p e is generating the induced macroscopic eddy currents, therefore, its distribution resembles the distribution of E (or j) in Figs. 11 i) and j). This instantaneous power is always positive, which implies that all the energy needed to generate these phenomena is, irreversibly, converted into heat. Furthermore, it is apparent that p e is highest at the surfaces and lowest in the center of the sheet. This corresponds to the distribution of j that is, theoretically, always zero in the center of the sheet. In contrast to this, the instantaneous power p h can be both positive as well as negative, which implies that the energy stored in the magnetic field inside the steel sheet flows partially back to the source. However, the average power over the whole period is still positive, which corresponds to the power loss due to static hysteresis. The results in Fig. 11 m) show that the highest power p h is required near the surfaces of the sheets, because the magnetic field is forced to these regions by the induced eddy currents. With evolving time in each half period, more power is transfered to the center of the sheet, which corresponds to the penetration of a decaying electromagnetic wave from the surfaces to the center of the sheet.
The viscous instantaneous power p v is, like p e , always positive and represents instantaneous power loss due to microscopic phenomena. According to (39), p v depends on dB s dt and H vs inside individual slices, therefore, the profiles shown in Figs. 11 q) and r) strongly resemble Figs. 11 e) and f). This power loss is also the highest near the surfaces of the steel sheets, because, in these regions, the rates of change of magnetic flux densities are the highest.
Finally, the instantaneous total power p t is shown in Figs. 11 s) and t). This power is the sum of all three discussed components, where the average value over a time period represents the total power loss inside the sheet. In the presented case, due to skin effect, substantially more power is dissipated near the edges of the sheet. Because pulsating magnetization is analyzed, the power loss is also pulsating, and is highest in the middle of each half period, where the highest rates of change of magnetic flux density inside the steel sheet occur.

2) SATURATED ELECTRICAL STEEL (B max = 1.5 T)
The influence of saturation is presented in Fig. 12. In this case the macroscopic eddy currents still generate a distinct phase shift betweenB s inside individual slices when the sheets are not saturated, as shown in Fig. 12 a). The magnetization begins at the surface, whereas a wave penetrates the sheet, as shown in Fig. 12 b). The penetrating wavefront has a very steep, almost rectangular shape, which corresponds to the non-linear magnetization property of the steel sheets.B s in individual slice changes quickly when the observed slice is not saturated, but is finally limited by saturation, as shown in Figs. 12 a) and e). When the steel sheet saturates, the electromagnetic field inside the sheet becomes almost spatially-uniformly distributed, where B max ≈ B max1 ≈ B max10 . When the steel sheet is saturated, H s inside the sheet starts to increase very rapidly, as shown in Figs. 12 c) and d). It is important to note that the dynamic effects are always influencing the distribution of the magnetic filed inside the steel sheets. However, in saturation, the penetration depth δ increases significantly due to low permeability, therefore, the skin effect can be neglected. Consequently, the spatial magnetization profiles of B(x/b) and H (x/b) are at material saturation almost constant functions, as shown in Figs. 12 b) and d).
The gradual spatial magnetization of the sheet can also be observed clearly from the rates of changes dB s dt that are shown in Figs. 12 e) and f). In contrast to the unsaturated case, dB s dt inside the individual slices have comparable maximum values. However, they are also distinctively shifted in phase. Furthermore, such a magnetization translates into gradual spatial generation of the viscous magnetic field H v , as shown in Figs. 12 g) and h). More interesting consequences of saturation are observed by analyzing the induced electric field E shown in Figs. 12 i) and j), as well as the magnetic field H e shown in Figs. 12 k) and l). Similar to the unsaturated case, the induced E is distributed non-linearly over the sheet's thickness when the sheet is not saturated. However, when all slices in the sheet are saturated, E is increasing linearly from the center to the edges of the sheet. This corresponds to the increased penetration depth and no skin effect, as discussed in the previous paragraph. When analyzing H e , shown in Figs. 12 k) and l), it can be observed that the impact of dB s dt in the inner slices on H e is significantly higher compared to the impact of dB s dt in the outer ones. This corresponds to (29) or (31), where it is apparent that the influence of dB s dt is multiplied in inner slices. Therefore, H e in Fig. 12 k) is increasing from the start of the half period and is the highest at t/T p ≈ 0.35, where the dB s dt are the highest in inner slices at the center of the sheet. At this time, dB s dt in outer slices are already significantly lower. This effect causes that the dynamic hysteresis loops inflate gradually with ongoing magnetization (in the presented case as they approach saturation). Therefore, the inclination of the dynamic loops is decreased with increasing excitation dynamics, as shown in Fig. 7 or Fig. 8. It is important to note that this effect is also visible in the unsaturated case presented in Figs. 11 e) and k), where substantially lower dB s dt in inner slices influence the generated H e substantially more compared to higher values dB s dt in outer slices.
Saturation has a significant influence on the power-loss distribution across the steel sheets. When all the slices in the sheet are saturated, the dynamic loss components p e and p v are the lowest, as presented in Figs. 11 m), n), q) and r). However, when the slices are de-saturated due to magnetization in the opposing direction, these loss components increase significantly. p e is still higher near the surfaces of the sheet, but penetrates slightly deeper with ongoing magnetization. However, in the center of the sheet, p e is still zero. In contrast to this, p v reflects the gradual magnetization of the sheet. Therefore, this loss component is at the beginning generated mainly near the surfaces, but, with ongoing magnetization, gradually penetrates the sheet, where the amplitude remains approximately constant. Right before the sheet is saturated again, this penetration reaches the center of the sheet and eventually decreases when the sheet is fully saturated. Similar to p v , p h is also reflecting the gradual spatial magnetization of the sheet, as presented in Figs. 11 o) and p). At the beginning of each half period the sheet is de-saturated, where a part of the stored energy in the magnetic field is returning to the source. Later in the magnetization process, however, energy is again required to magnetize the sheet in the opposite direction. Because of the gradual magnetization process, the peak value of p h travels gradually to the center of the sheet, whereas increased power p h is needed because the sheet is saturated in the opposite direction. Because the static magnetization process consists of irreversible and reversible magnetization processes, a part of the supplied energy is stored in the magnetic field, and, in the next magnetization cycle, is returned to the power source. Static hysteresis power loss can be obtained by calculating the average value of p h over one magnetization cycle. The time and space distributions of p h are presented in Figs. 11 s), and t). The obtained results show that the gradual magnetization process, in the beginning, requires more power at the surfaces of the sheet, whereas the power is later transferred into the sheet's thickness. Compared to the unsaturated case, the power, as well as the power loss, is distributed more uniformly across the sheet's thickness.

V. CONCLUSION
This paper starts with a thorough introduction to the fundamentals of non-linear eddy-current diffusion in NO soft magnetic materials. The importance is pointed out of the specific material structure of NO electrical steels that enables straightforward application of Maxwell equations. Based on the specific geometry of steel sheets and proper boundary conditions, a 1-D lamination model with hysteresis was derived and transferred to the theory of the PMD model. The basic (classical) PMD model, which accounts for hysteresis and non-local (macroscopic) eddy-current effects, was then extended to consider the effect of local (microscopic) eddy currents by means of magnetic viscosity. Both versions of the model were parameterized by using standardized measured data, and compared to measured data at high frequencies under sinusoidal magnetization waveforms. Two distinctively different NO electrical steel grades, i.e., in terms of iron thickness, silicon content, grain size, were studied. A very good accuracy of the model, despite the application to different materials and merely physical-based parameters, was obtained. The importance of considering non-linear, hysteretic properties and the magnetic viscosity when calculating dynamic variables and power loss inside NO soft magnetic steel sheets with non-uniform magnetic fields is pointed out. Both PMD model versions were compared directly in terms of identification procedure facilities, accuracy, spatial discretization and numerical implementation.
The application of the coupled approach (lamination model plus static hysteresis model plus magnetic viscosity) to unidirectional, pulsating high frequency excitation cases provides a detailed insight and a thorough analysis of the temporal and spatial distribution of power-loss and field components in the electrical steel sheets. Further on, the separation of the iron-loss in the individual components allows one to disentangle the closely intertwined effects at different spatial scales and, ultimately, to draw a link to the results obtained from statistical theory of losses. The ability to predict the dynamic hysteresis loop shape, electromagnetic variables and power-loss components is not limited to harmonic excitation cases, but also enables the analysis of arbitrary distorted waveforms. The presented results enabled an application-specific selection of the lamination model, taking account of the actual microstructural and geometric properties of the soft magnetic material. Therefore, error-prone and limited engineering approaches commonly used for iron-loss calculation based on a huge amount of measured data can be replaced by the PMD model, coupled with the most suited hysteresis model for the considered operation conditions. As a result, not only the hysteresis loop shape and magnetization dynamics can be calculated, but also the parameter identification effort is reduced to two quasi-static (low-frequency) measurements or one static to identify the major hysteresis loop and the usage of material specific data, e.g., specific electrical conductivity and thickness.
Based on the presented analysis, it was shown that the potentially biggest improvement of the discussed models could be achieved with further research and development of adequate static hysteresis descriptions that can describe arbitrary magnetization profiles. Future research work will, therefore, be focused on the development of efficient and accurate vector hysteresis models that will enable modeling of the anisotropic magnetic properties of NO materials. This will, in turn, enable adequate extension of the discussed 1-D lamination model to the 2-D excitation cases, i.e., accurate modeling of magnetization dynamics and power loss of rotating magnetic fields.

APPENDIX DETAILED CALCULATION OF THE PMD MODEL A. DISTRIBUTION OF THE INDUCED ELECTRIC FIELD
Function E s (x), described by (21), can be interpreted alternatively as the sum of changes of the induced electric field E inside the slices by In (45), to the known BC E (x = 0) = 0, all total changes of induced electric field E s in all the inner slices are added and, finally, the position x dependent change of induced electric field E s (x) in the slice s is added. The latter is determined by where E s (x) is changing linearly across each slice. By substituting x with the upper border position of the observed slice x = x us , the total change of induced electric field E s across arbitrary slice s is determined by Distribution function E s (x), described by (22), is obtained by applying (46) and (47) to (45).

B. DISTRIBUTION OF THE MAGNETIC FIELD
Analogous to (45), H s (x), described by (26), can be interpreted as the sum of known Dirichlet BC at the surface (i.e. H sur ), changes of magnetic field H s across all the outer slices, and change of the magnetic field H s (x) from the upper border x us to arbitrary position x inside the slice s. This is described by H s (x) can be determined by calculating the third term on the right hand side in (26), where the approximated distribution of eddy-current density j s (x) described by (24) is applied. Because j s (x) is already represented as a sum, H s (x) is obtained by integrating all its terms by The result of integration in (49) is given by By substituting x in (50) with the upper border position of the observed slice x = x us , the total change of the magnetic field H s across arbitrary slice s is determined by The obtained expressions (50) and (51) can be applied to (48), whereas H s (x) can be expressed by The second term on the right hand side of (52) is described by a sum of the sums (i.e. a double sum), which can be expressed as two single sum terms by H s (x) expressed by (27) is obtained by applying (53) in (52).