Approximation of Hysteresis Changes in Electrical Steel Sheets

: This paper describes a simple method of approximating hysteresis changes in electrical steel sheets. This method is based on assumptions that ﬂux density or ﬁeld strength changes are a sum or a difference of functions that describe one curve of the limiting hysteresis loop and a certain ‘transient’ component. Appropriate formulas that present the ﬂux density as functions of the ﬁeld strength and those that present inverse dependencies are proposed. An application of this approximation requires knowledge of the measured limiting hysteresis loop and a few minor loops. Algorithms for determining changes in the ﬂux density or ﬁeld strength are proposed and discussed. The correctness of the proposed approximation of hysteresis changes was veriﬁed through a comparison of measured hysteresis loops with the loops calculated for several different excitations of the magnetic ﬁeld occurring in dynamo and transformer steel sheets. Additionally, an example of the application of the proposed approximation of hysteresis changes is discussed in the paper. The proposed approximation of hysteresis changes is recommended for numerical calculations of the magnetic ﬁeld distribution in dynamo and transformer steel sheets.


Introduction
The hysteresis phenomenon is a characteristic feature of any magnetic circuit containing a ferromagnetic material. This phenomenon is related to a certain delay in changes in the magnetic flux density in relation to changes in the magnetic field strength in this material. The hysteresis mechanism is quite well known [1][2][3][4][5][6][7]; however, the formulation of an appropriate mathematical model based on the knowledge of the phenomena describing the magnetization process is still a complex issue. Therefore, for many decades, intensive research has been carried out to develop a hysteresis model that would be sufficiently accurate and useful in the field calculations of ferromagnetic materials.
On a micro scale, the study of magnetization processes is often based on the Landau-Lifschitz-Gilbert equation, which describes the behavior of a single magnetic dipole in an external magnetic field [4][5][6]8]. The original equation formulated by Landau-Lifschitz was modified by Gilbert to take into account the damping conditions of the magnetic moment rotation, which in the steady state will occupy a position along the direction of the magnetic field strength. As stated in a previous study [9], the calculation results obtained on the basis of the Landau-Lifschitz-Gilbert equation "can be averaged to provide hysteresis curves of material". However, it should be mentioned that hysteresis curves of ferromagnetic materials are also often described by different empirical fitting procedures.
The best known models, which are based on energy relationships, are the Stoner-Wohlfarth and Jiles-Atherton models. The former presents the magnetization process as rotations of magnetic moments of non-interacting single-domain particles [10][11][12]. To take into account the interaction of particles, the concept of the effective field has been introduced, which is the sum of the field strength of the external magnetic field and the product of the magnetization and a certain dimensionless mean field parameter [4,8,[10][11][12][13].
The resultant magnetization vector in the Stoner-Wohlfarth model is a sum of the magnetizations of single particles aligned along specified directions. However, displacements of the domain wall are not considered. The application of this model for macroscopic samples must specify a certain number of directions for a particular ferromagnetic material. Therefore, computation times are relatively long, and they depend directly on the number of specified directions.
The assumption that the energy supplied to a sample of an isotropic polycrystalline material is a sum of the change in the magnetostatic energy and hysteresis loss is the basis for the Jiles-Atherton model [5,9,11,12]. This model requires five physical parameters to describe hysteresis curves. However, in many instances, these parameters must be modified to obtain an effective representation of the measured waveforms. One difficulty is the determination of an anhysteretic curve because this curve is a theoretical curve. The Jiles-Atherton model is used in modelling relatively narrow hysteresis loops. Frequently, the model parameters should be determined for almost each amplitude of the field strength. This causes a certain inconvenience when using this model; nevertheless, this model and its modifications are still used [14][15][16][17][18][19][20][21][22][23][24][25].
The phenomenological Preisach model assumes that a considered ferromagnetic material consists of particles (hysterons), which are magnetized to a positive or negative saturation. A characteristic feature of this model is the Preisach triangle consisting of the following two parts: the first part refers to positively magnetized particles, and the second part refers to negatively magnetized particles. A certain distribution function, whose determination is time-consuming, describes how many hysterons are located at particular points of the Preisach triangle. This model enables the creation of any hysteresis or minor loop, depending on an appropriate distribution function. However, long calculation times and dependence on the degree of Preisach triangle discretization are significant disadvantages of this model. The Preisach model is still modified and often used in different applications [26][27][28][29][30][31][32][33][34][35][36].
Hysteresis models formulated with the use of special operators, the so-called 'play' and 'stop' hysterons, are also phenomenological [37]. These types of models were developed by Matsuo et al. [38][39][40][41]. Some hysteresis models are based on certain differential equations that relate the magnetic field strength to the magnetic flux density. However, these equations do not result directly from the mechanism of phenomena occurring in the magnetization process. The Hodgdon model belongs to this category [11,42]. The determination of functions occurring in the Hodgdon model is time-consuming and requires a series of tests. Note that slight changes in the parameters of these functions significantly affect the shape of the calculated hysteresis loops.
There are other, less known hysteresis models in the literature. One of them is the Globus model, which assumes that the behavior of the entire ferromagnetic sample can be reduced to one spherical grain [12]. This group also includes the hysteresis model developed by Chua [11]. In this model, the resultant magnetic field strength is the sum of the field strength resulting from the hysteresis-free magnetization curve and the field strength dependent on the derivative of the magnetic flux density versus time; in fact, the Chua model is a dynamic model. To describe changes in electrical steel sheets, the Pry and Bean domain model is also applied [43,44].
Calculations of the magnetic field distribution in electrical steel sheets often require consideration of the hysteresis phenomenon. These sheets are used in the construction of the magnetic cores of transformers, rotating machines, and other electrical devices. In many scenarios, a certain model of the hysteresis phenomenon must be considered when calculating the magnetic field distribution and hysteresis loss of a tested magnetic circuit. Instead of one unambiguous relationship, these models propose different dependencies between the flux density and field strength as the hysteresis changes. However, the formulation of a suitable mathematical hysteresis model based on the hysteresis mechanism is still a relatively difficult problem. Some magnetic hysteresis models refer to the so-called dynamic hysteresis loop, which considers the effects of eddy currents on changes in Energies 2021, 14, 4110 3 of 18 flux density. However, in many proposals, eddy currents are considered using separate differential equations [45,46].
It should be clearly emphasized that the purpose of this paper is not to formulate a new hysteresis model but rather to present a proposal for approximation of hysteresis changes; this approximation should be easily applicable to equations describing the magnetic field distribution. In calculations, the area of the magnetic field is divided into several dozen or even more thousand elementary segments. Owing to the nonlinear nature of electrical sheets, the magnetic flux density in individual segments changes in different ways. Therefore, the hysteresis phenomenon should be included in all elementary segments concerning a considered ferromagnetic material.

Approximation of Flux Density Changes
Each hysteresis model briefly discussed in the previous chapter has advantages and disadvantages. However, note that any model of the magnetic hysteresis should have the following properties: • Irreversible displacements of domain walls are the main cause of the hysteresis phenomenon in soft magnetic materials; these displacements are caused by changes in field strength in a sample of ferromagnetic materials. When the field strength values are low, then the magnetization process is reversible. However, for most soft ferromagnetic materials, particularly for electrical steel sheets, the effect of the reversible magnetization can be neglected. Notably, in a high-value range of the field strength, the magnetization process is also reversible, but then this process refers to reversible rotations of magnetization vectors towards the external magnetic field.
It follows from the nature of the hysteresis phenomenon that any point P with coordinates (H 0 , B 0 ) inside the limiting hysteresis loop moves along a certain trajectory to the bottom or upper limiting curve depending on changes in the field strength ( Figure 1). This means that the difference B 0 − B b (H 0 ) between the initial flux density B 0 at point P and the flux density value B b (H 0 ) on the bottom limiting curve determined for H 0 should decrease to zero if the field strength H increases. Therefore, the difference B 0 − B b (H 0 ) can be treated as an 'initial value' of the certain transient component of changes in the flux density of point P moving to the bottom limiting curve B b (H).
In this paper, the subscripts B and H concern parameters which are related to changes in the flux density or magnetic field, respectively; subscripts b and u refer to the bottom or upper curve of the hysteresis loop, respectively, and the subscript r is used if the magnetic flux or field strength increases, otherwise the subscript d is used; for example, the parameter a Bur relates to the case when the flux density B is a function of the field strength, u refers to the upper curve of the limiting hysteresis loop, and r is used when the flux density increases; in turn, the parameter a Hbd relates to the case when the field strength H depends on the flux density, b refers to the bottom curve of the limiting hysteresis loop, and d is used when the field strength decreases. where Bb(H) is the bottom curve of the limiting hysteresis loop, ΔBr = B0 − Bb(H0), H0 and B0 denote the coordinates of point P(H0,B0), and kBr is an attenuation coefficient of the transient component when H increases.
The second component of Equation (1) is a certain pseudotransient component of the changes in flux density. The coefficient kBr(H) affects the 'dynamics' at which point P moves to the bottom limiting hysteresis curve. Generally, this coefficient is a nonlinear function of H. When H decreases, the difference Bu(H0) − B0 between the flux density value Bu(H0) on the upper limiting curve, determined for H0, and the initial flux density B0 at point P should decrease to zero. Therefore, the difference Bu(H0) − B0 can be treated, similarly to previously, as a certain initial value of the transient component of changes in the flux density of point P moving to the upper limiting curve Bu(H). In this case, the change in the flux density Bd(H) of point P can be approximated as follows: where To apply this description of hysteresis changes, the limiting hysteresis loop must be measured for a considered electrical sheet. Here, when only the coercive field, residual flux density, and saturation flux density are known, curves of the limiting hysteresis loop can be approximated using, for example, an arctg function which considers the abovementioned parameters.
The coefficients kBr and kBd depend on the position of point P(H0,B0). Owing to the nature of the hysteresis phenomenon, the derivative dBr/dH on the trajectory of point P(H0,B0) for increasing values of H is in the range (aBur,aBbr), where aBur describes the value of this derivative when point P(H0,B0) is situated on the upper limiting curve Bu(H) ( Figure  1), and aBbr denotes this derivative when point P(H0,B0) is situated sufficiently close to the bottom hysteresis curve, i.e., to the curve to which it moves when H increases. Thus, aBbr is equal to the derivative dBb/dH of the bottom curve at point Pr ( Figure 1). The value aBur is determined based on several measured symmetrical hysteresis loops. Denoting the derivative dBr/dH, whose value depends on the position of the point P(H0,B0) as aBr, it can be expressed as follows: Using some analogy to transient states occurring in electrical circuits, where the transient component of a certain quantity decreases and its changes are described by means of an exponential function, the flux density changes of point P(H 0 ,B 0 ) can be described as follows when the field strength H increases: When H decreases, the difference B u (H 0 ) − B 0 between the flux density value B u (H 0 ) on the upper limiting curve, determined for H 0 , and the initial flux density B 0 at point P should decrease to zero. Therefore, the difference B u (H 0 ) − B 0 can be treated, similarly to previously, as a certain initial value of the transient component of changes in the flux density of point P moving to the upper limiting curve B u (H). In this case, the change in the flux density B d (H) of point P can be approximated as follows: To apply this description of hysteresis changes, the limiting hysteresis loop must be measured for a considered electrical sheet. Here, when only the coercive field, residual flux density, and saturation flux density are known, curves of the limiting hysteresis loop can be approximated using, for example, an arctg function which considers the abovementioned parameters.
The coefficients k Br and k Bd depend on the position of point P(H 0 ,B 0 ). Owing to the nature of the hysteresis phenomenon, the derivative dB r /dH on the trajectory of point P(H 0 ,B 0 ) for increasing values of H is in the range (a Bur ,a Bbr ), where a Bur describes the value of this derivative when point P(H 0 ,B 0 ) is situated on the upper limiting curve B u (H) (Figure 1), and a Bbr denotes this derivative when point P(H 0 ,B 0 ) is situated sufficiently close to the bottom hysteresis curve, i.e., to the curve to which it moves when H increases. Thus, a Bbr is equal to the derivative dB b /dH of the bottom curve at point P r (Figure 1). The value a Bur is determined based on several measured symmetrical hysteresis loops. Denoting the derivative dB r /dH, whose value depends on the position of the point P(H 0 ,B 0 ) as a Br , it can be expressed as follows: where Figure 1). If the field strength H increases and the initial position of point P is close to the upper limiting curve B u (H), then a Br is almost equal to a Bur because d Br is close to 0. In turn, when point P is close to the bottom curve B b (H), d Br tends to 1, and a Br is almost equal to a Bbr . This means that point P moves along the bottom limiting curve B b (H); these are extreme cases resulting from Equation (3).
The dimensionless exponent p B affects changes in a Br , and it should be selected in such a manner that the differences between the calculation results and real transients are the lowest; frequently, this exponent is higher than 1. The exponent p B has higher values when the shape of the hysteresis loop is similar to a rectangle. Note that the exponent p B should have the same value independently on the field strength changes.
To calculate the flux density value B r (H) for the next increasing field strength H, it is necessary to determine the value of the coefficient k Br based on the flux density B obtained from the previous solution. Assuming that the coefficient k Br has a constant value in a sufficiently small neighborhood of the value H 0 , then, based on Equation (1), the following relationship can be expressed for H = H 0 : Because the derivative dB r (H 0 )/dH, determined for H 0 , is equal to a Br , the coefficient k Br (H) is determined using the following expression: The coefficient a Br for point P(H 0 ,B 0 ) is calculated based on Equation (3), and then the coefficient k Br (H) is determined using Equation (5). The value of this coefficient is inserted into Equation (1), and for the next H, the new flux density B r (H) is calculated.
Denoting the derivative dB d /dH, whose value also depends on the position of the point P(H 0 ,B 0 ), as a Bd (similar to a Br ), it can be written for decreasing field strength as follows: where , and a Bud is equal to derivative dB d /dH when the point P(H 0 ,B 0 ) is located sufficiently close to the upper hysteresis curve (point P d in Figure 1). For example, Figure 2 shows the dependencies of the parameters a Br and a Bd on H for B 0 = 0 determined for dynamo sheet M530-50A and transformer sheet M120-27S; for this value of B 0 , the field strength can change in the range (−H c ,H c ). Note that the curves shown in Figure 2 are similar; however, the ranges of the field strength and flux density changes are different.
With such a simplification, it is only necessary to determine the value of the parameter pB and to introduce the functions describing the limiting hysteresis loop of a considered electrical steel sheet into the calculation algorithm. It is advantageous when it can be assumed that the coefficients kBr and kBd have constant values for any external excitations; then, the calculation times are shorter, but differences between the measured and approximated flux density values are greater than in the case when these coefficients change their value depending on H. The distances between point P and the limiting hysteresis curves are determined for initial point P(B0,H0). Depending on changes in H, aBr and aBd are calculated. Next, the coefficients kBr or kBd are determined, and a new value of the flux density Br or Bd is calculated. The algorithm for determining changes in flux density is shown in Figure 3.
The curves of the limiting hysteresis loop (described by appropriate analytical functions or presented by means of array functions) of the considered steel sheet should be introduced into the algorithm. The coefficient aBur (aBbd is equal to aBur) is determined on the basis of several measured minor loops, or these coefficients are assumed to be equal to zero. The parameter pB is selected by comparing several measured loops with approximated curves. Depending on changes in the field strength, the flux density value Br(H) or Similarly, to calculate the flux density value B d (H), according to Equation (2), for the next decreasing field strength H, it is necessary to determine the coefficient k Bd based on the flux density B obtained from the previous solution. Assuming that the coefficient k Bd has a constant value in a sufficiently small neighborhood of the value H 0 , then, based on Equation (2), the following relationship can be written for H = H 0 : The derivative dB d (H 0 )/dH, determined for H 0 , is equal to a Bd , so the coefficient k Bd (H) is determined as follows: The coefficient a Bd for point P(H 0 ,B 0 ) is calculated based on Equation (6), and then the coefficient k Bd (H) is determined using Equation (8). The value of this coefficient is inserted into Equation (2), and for the next decreasing value H, the new flux density B d (H) is calculated.
Note that the coefficients k Br and k Bd are not selected a priori, but they are numerically determined for the values a Br and ∆B r or a Bd and ∆B d determined in the previous calculation step. Before calculations, the value of the exponent p B must be selected; the main criterion is the comparison between calculated minor loops for several ranges of the field strength changes and the corresponding measured loops.
It is worth noting that the parameter a Bbd is equal to parameter a Bur ; it significantly simplifies the procedure of parameter determination of the proposed approximation of hysteresis changes. Numerical calculations performed for different electrical sheets and different excitations indicate that satisfactory results are obtained assuming that parameters a Bur and a Bbd are equal to zero; then, Equations (3) and (6) simplify to the following form: With such a simplification, it is only necessary to determine the value of the parameter p B and to introduce the functions describing the limiting hysteresis loop of a considered electrical steel sheet into the calculation algorithm. It is advantageous when it can be assumed that the coefficients k Br and k Bd have constant values for any external excitations; then, the calculation times are shorter, but differences between the measured and approximated flux density values are greater than in the case when these coefficients change their value depending on H.
The distances between point P and the limiting hysteresis curves are determined for initial point P(B 0 ,H 0 ). Depending on changes in H, a Br and a Bd are calculated. Next, Bd(H) is determined for the new value of H. Note that the proposed approximation requires determining the limiting hysteresis curves and only two parameters aBur and pB, unlike, e.g., the Jiles-Atherton model for which five parameters must be determined.

Approximation of Field Strength Changes
In numerical calculations particularly involving magnetic field distributions in electrical steel sheets, an inverse approximation H = f(B) of the hysteresis loop is frequently applied. In those calculations, the vector potential is an unknown variable in the equation system describing the magnetic field distribution. In this approximation, the bottom lim- When B increases, the equation describing field strength changes can have the following form: where  The curves of the limiting hysteresis loop (described by appropriate analytical functions or presented by means of array functions) of the considered steel sheet should be introduced into the algorithm. The coefficient a Bur (a Bbd is equal to a Bur ) is determined on the basis of several measured minor loops, or these coefficients are assumed to be equal to zero. The parameter p B is selected by comparing several measured loops with approximated curves. Depending on changes in the field strength, the flux density value B r (H) or B d (H) is determined for the new value of H. Note that the proposed approximation requires determining the limiting hysteresis curves and only two parameters a Bur and p B , unlike, e.g., the Jiles-Atherton model for which five parameters must be determined.

Approximation of Field Strength Changes
In numerical calculations particularly involving magnetic field distributions in electrical steel sheets, an inverse approximation H = f (B) of the hysteresis loop is frequently applied. In those calculations, the vector potential is an unknown variable in the equation system describing the magnetic field distribution. In this approximation, the bottom limiting hysteresis curve in the approximation B = f (H) becomes the upper limiting curve where Hb(B) is the function describing the lower curve of the limiting hysteresis loop, ΔHd = H0 − Hb(B0), and kHd can be considered an attenuation coefficient of a 'transient' component when B decreases. Assuming that aHr is equal to the derivative dHr/dB, the coefficient aHr can be expressed as follows: where dHr = dH1/(dH1 + dH2) (Figure 4). Note that in this approach, H0 − Hb(B0) is identically dH1, and Hu(B0) − H0 is identically dH2.
If B increases and the initial position of point P is close to the bottom limiting curve Hb(B), then aHr is almost equal to aHbr because dHr is close to 0. In turn, when point P is close to the upper limiting curve Hu(B), then dHr tends to 1 and aHr is almost equal to aHur; then, point P moves along the upper curve Hu(B). The exponent pH in Equation (12) is less than 1, and it should be experimentally determined for a tested electrical steel sheet. To calculate the field strength Hr(B) for the next increasing flux density B, it is necessary to determine the value of the coefficient kHr based on the value H obtained from the previous solution.
Assuming that kHr has a constant value in a certain sufficiently small neighborhood of the flux density B0, then based on Equation (10), the following relationship can be obtained: The derivative dHr(B0)/dB is equal to aHr; thus, the coefficient kHr can be calculated as follows: When B increases, the equation describing field strength changes can have the following form: where H b (B) is the function describing the lower curve of the limiting hysteresis loop, Assuming that a Hr is equal to the derivative dH r /dB, the coefficient a Hr can be expressed as follows: dH r dB = a Hr = a Hbr + (a Hur − a Hbr ) (d Hr ) p H , where d Hr = d H1 /(d H1 + d H2 ) (Figure 4). Note that in this approach, H 0 − H b (B 0 ) is identically d H1 , and H u (B 0 ) − H 0 is identically d H2 .
If B increases and the initial position of point P is close to the bottom limiting curve H b (B), then a Hr is almost equal to a Hbr because d Hr is close to 0. In turn, when point P is close to the upper limiting curve H u (B), then d Hr tends to 1 and a Hr is almost equal to a Hur ; then, point P moves along the upper curve H u (B). The exponent p H in Equation (12) is less than 1, and it should be experimentally determined for a tested electrical steel sheet. To calculate the field strength H r (B) for the next increasing flux density B, it is necessary to determine the value of the coefficient k Hr based on the value H obtained from the previous solution.
Assuming that k Hr has a constant value in a certain sufficiently small neighborhood of the flux density B 0 , then based on Equation (10), the following relationship can be obtained: The derivative dH r (B 0 )/dB is equal to a Hr ; thus, the coefficient k Hr can be calculated as follows: The value a Hr for the point P(H 0 ,B 0 ) is determined based on Equation (12), and the next value of the coefficient k Hr is calculated using Equation (14). This value is considered in Equation (10), and a new field strength value H r (B) is calculated for the next value B.
When B decreases, the coefficient a Hd can be expressed as follows: where d Hd = d H2 /(d H1 + d H2 ) and a Hbd is equal to the derivative dH d /dB when point P(H 0 ,B 0 ) is located on the bottom limiting curve H b (B).
To calculate the field strength H d (B) for the next decreasing flux density B, the coefficient k Hd should be determined based on the field strength H obtained from the previous solution. Assuming that the coefficient k Hd has a constant value in a sufficiently small neighborhood of the value B 0 , then, based on Equation (11), the following relationship can be written for B = B 0 : In this case, the derivative dH d (B 0 )/dB is equal to a Hd , so the coefficient k Hd can be determined as follows: In this approximation, the parameters a Hbr and a Hud are also equal to each other (similar to the previous approximation); however, they cannot be omitted as in the approximation B = f (H). In practice, they are equal to the reciprocal of the nonzero parameter a Bur occurring in Equation (3).
The distances between point P and the limiting hysteresis curves are calculated for the given position of point P(B 0 ,H 0 ). The coefficients a Hbr and a Hud are determined depending on the flux density changes, and a Hr and k Hr or a Hd and k Hd are calculated. This enables the calculation of a new field strength value H r or H d . The dependencies of parameters a Hr and a Hd on B for H 0 = 0 are shown in Figure 5. Analogous to Figure 2, the curves shown in Figure 5 are similar in shape; however, the ranges of changes in both the flux density and field strength are different. The algorithm to determine the field strength changes is similar to that shown in Figure 6. the given position of point P (B0,H0). The coefficients aHbr and aHud are determined depending on the flux density changes, and aHr and kHr or aHd and kHd are calculated. This enables the calculation of a new field strength value Hr or Hd. The dependencies of parameters aHr and aHd on B for H0 = 0 are shown in Figure 5. Analogous to Figure 2, the curves shown in Figure 5 are similar in shape; however, the ranges of changes in both the flux density and field strength are different. The algorithm to determine the field strength changes is similar to that shown in Figure 6.      Figure 7 shows a family of first-reversal curves calculated for increasing and decreasing values of the field strength in transformer sheet M120-27S (magnetization along the rolling direction); qualitatively similar curves have the non-oriented dynamo sheet M530-50A, the difference mainly concerns coercivity values. Regardless of the initial position of point P, the trajectories of points tend to the bottom or upper limiting curve; this conclusion refers to both approximations.   Figure 9 shows examples of minor loops when H has a constant component. Regardless of the initial point P1 or P2, the minor loops should be the same in steady states. The minor loops reach a final course after several cycles of changes in the field strength. Note that the minor loops can be determined using both the Preisach and Jiles-Atherton models; however, the latter model must be extended to calculate minor hysteresis loops. Each trajectory of point P (Figures 1 and 4) is uniquely defined by the initial location of point P; thus, the first Madelung's rule is considered to be satisfied [47,48]. The second Madelung's rule states that each hysteresis loop should be closed (so-called 'return-pointmemory'). Corresponding minor hysteresis loops must be symmetrical with respect to the coordinate system center. However, this scenario occurs only in the steady state of hysteresis changes (Figure 8) because minor loops drift gradually to these closed loops [30,31]. The closed minor loops calculated for the assumed field strength or flux density changes should be the same independent of the starting point ( Figure 9); this means that the third Madelung's rule, the so-called congruency property, is fulfilled. Additionally, the fourth rule is satisfied because the hysteresis loops are symmetrical with respect to the coordinate system center for alternating changes in the flux density.  Figure 9 shows examples of minor loops when H has a constant component. Regardless of the initial point P1 or P2, the minor loops should be the same in steady states. The minor loops reach a final course after several cycles of changes in the field strength. Note that the minor loops can be determined using both the Preisach and Jiles-Atherton models; however, the latter model must be extended to calculate minor hysteresis loops. Each trajectory of point P (Figures 1 and 4) is uniquely defined by the initial location of point P; thus, the first Madelung's rule is considered to be satisfied [47,48]. The second Madelung's rule states that each hysteresis loop should be closed (so-called 'return-point-memory'). Corresponding minor hysteresis loops must be symmetrical with respect to the coordinate system center. However, this scenario occurs only in the steady state of hysteresis changes ( Figure 8) because minor loops drift gradually to these closed loops [30,31]. The closed minor loops calculated for the assumed field strength or flux density changes should be the same independent of the starting point ( Figure 9); this means that the third Madelung's rule, the so-called congruency property, is fulfilled. Additionally, the fourth rule is satisfied because the hysteresis loops are symmetrical with respect to the coordinate system center for alternating changes in the flux density.

Properties of the Proposed Approximation Methods
trajectory of point P (Figures 1 and 4) is uniquely defined by the initial location of point P; thus, the first Madelung's rule is considered to be satisfied [47,48]. The second Madelung's rule states that each hysteresis loop should be closed (so-called 'return-pointmemory'). Corresponding minor hysteresis loops must be symmetrical with respect to the coordinate system center. However, this scenario occurs only in the steady state of hysteresis changes (Figure 8) because minor loops drift gradually to these closed loops [30,31]. The closed minor loops calculated for the assumed field strength or flux density changes should be the same independent of the starting point ( Figure 9); this means that the third Madelung's rule, the so-called congruency property, is fulfilled. Additionally, the fourth rule is satisfied because the hysteresis loops are symmetrical with respect to the coordinate system center for alternating changes in the flux density.

Approximations of Hysteresis Changes
Dynamo and transformer steel sheets are typical ferromagnetic material characterized by the hysteresis loop. Depending on assumed unknowns occurring in equations of the magnetic field distribution, the basic hysteresis model or its inverse version is applied. However, note that the calculation results of the hysteresis changes should be the same using both versions of the hysteresis model. Numerical calculations of hysteresis changes were performed for the non-oriented dynamo sheet M530-50A and grain-oriented transformer sheet M120-27S.
The dynamo sheets M530-50A are produced as non-oriented steel sheets; however, most of these sheets have some anisotropic properties [49]. In turn, transformer sheets M120-27S are typical grain-oriented sheets. Owing to the occurrence of the Goss texture, these sheets have different magnetic properties in individual directions on the sheet plane. Therefore, a different hysteresis loop should be considered for each magnetization direction. This is significant in calculations of the magnetic field distribution at the corners of transformer cores and in areas where the columns connect to the transformer yokes.
It was assumed that the parameters a Bur and a Bbd concerning both considered sheets were equal to zero in the approximation B = f (H). However, in the inverse approximation H = f (B), the parameters a Hbr and a Hud were equal to 3000 A/Tm and 1500 A/Tm for the dynamo and transformer sheet, respectively. The hysteresis changes were approximated for several assumed maximum values H max of the field strength, selecting appropriate values of the parameters p B and p H (Table 1). It is worth noting that the parameter p B in the approximation B = f (H) has the same value, regardless of the value H max , while the value of the parameter p H depends on the value H max , especially when the flux density is close to the saturation value.  Figure 10 shows symmetrical hysteresis loops calculated using the basic and inverse approximation. Additionally, for comparison, the measured hysteresis loops are also presented in Figure 10. (a) (b) Figure 11. Hysteresis loops of transformer sheet M120-27S: (a) for the 45° direction; (b) for the transverse direction; line markings as in Figure 10.
The key meaning in the proposed approximation of hysteresis changes is to determine the exponents pB or pH by comparing the calculated changes in the flux density or field strength with the corresponding measured hysteresis loops. The assumption of the exponents pB and pH other than those given in Table 1, has an influence on the shape of minor loops. Figure 12 shows the minor loops of the transformer sheet, determined for the exponents whose values differ by 20% with respect to the values given in Table 1. The hysteresis loops of transformer sheets have a 'classical' shape when the magnetization processes occur along the rolling direction or near that direction (Figure 10b). When the magnetization process takes place at an angle greater than 45 degrees relative to the rolling direction, the shapes of the hysteresis loops vary considerably with respect to the classical hysteresis loop. In the corners and the so-called T-joint points in magnetic circuits of the three-phase transformer, the magnetization processes run along directions other than the rolling direction. Figure 11 shows hysteresis loops of the mentioned transformer sheet measured and calculated for the 45 • direction and for the transverse direction. (a) (b) Figure 11. Hysteresis loops of transformer sheet M120-27S: (a) for the 45° direction; (b) for the transverse direction; line markings as in Figure 10.
The key meaning in the proposed approximation of hysteresis changes is to determine the exponents pB or pH by comparing the calculated changes in the flux density or field strength with the corresponding measured hysteresis loops. The assumption of the exponents pB and pH other than those given in Table 1, has an influence on the shape of minor loops. Figure 12 shows the minor loops of the transformer sheet, determined for the exponents whose values differ by 20% with respect to the values given in Table 1. The key meaning in the proposed approximation of hysteresis changes is to determine the exponents p B or p H by comparing the calculated changes in the flux density or field strength with the corresponding measured hysteresis loops. The assumption of the exponents p B and p H other than those given in Table 1, has an influence on the shape of minor loops. Figure 12 shows the minor loops of the transformer sheet, determined for the exponents whose values differ by 20% with respect to the values given in Table 1. The key meaning in the proposed approximation of hysteresis changes is to determine the exponents pB or pH by comparing the calculated changes in the flux density or field strength with the corresponding measured hysteresis loops. The assumption of the exponents pB and pH other than those given in Table 1, has an influence on the shape of minor loops. Figure 12 shows the minor loops of the transformer sheet, determined for the exponents whose values differ by 20% with respect to the values given in Table 1.    Table 1, red and blue solid lines-values of the coefficients higher or lower 20% according to values in Table 1, respectively.
Using the described approximation of the hysteresis changes, the flux density and field strength changes for other ferromagnetic materials can be calculated. For example, Figure 13 shows the hysteresis loops of Mn-Mg ferrite and ferromagnetic material ЮН14ДK24, which is applied in hysteresis motors.
Energies 2021, 14, x FOR PEER REVIEW 14 of 17 Using the described approximation of the hysteresis changes, the flux density and field strength changes for other ferromagnetic materials can be calculated. For example, Figure 13 shows the hysteresis loops of Mn-Mg ferrite and ferromagnetic material ЮН14ДK24, which is applied in hysteresis motors.

Application Notes
In modeling magnetic phenomena, the time consumption is a significant problem. On the one hand, this problem concerns determination of the model parameters; on the other hand, how long are calculation times for the assumed changes in the field strength or flux density. In the proposed approximation, it is necessary to know the limiting hysteresis loop and to assume the value of the exponent the pB or pH, which can be obtained after a few comparisons of the calculated minor hysteresis loops with the measured loops. It should be noted that the Preisach model requires the determination of the so-called grain distribution function, and in the Jiles-Atherton model five parameters must be assumed, whose values depend on changes in the field strength. The requirement to know the limiting hysteresis loop and a few minor loops is unquestionable because a comparison between calculation and measurements results should be the basic criterion for assessing the correctness and usefulness of the proposed approximation of the hysteresis changes.
Calculation times of a single hysteresis loop are comparable. On the one hand, a comparison between these times using different models is not justified, because errors between the calculation and measurement results depend, inter alia, on the assumed number of calculation steps. On the other hand, it is more important how long the calculation times are when a hysteresis model is included in a large-dimensional system of matrix algebraic equations describing the magnetic field distribution. Often the computation time is of secondary im-

Application Notes
In modeling magnetic phenomena, the time consumption is a significant problem. On the one hand, this problem concerns determination of the model parameters; on the other hand, how long are calculation times for the assumed changes in the field strength or flux density. In the proposed approximation, it is necessary to know the limiting hysteresis loop and to assume the value of the exponent the p B or p H , which can be obtained after a few comparisons of the calculated minor hysteresis loops with the measured loops. It should be noted that the Preisach model requires the determination of the so-called grain distribution function, and in the Jiles-Atherton model five parameters must be assumed, whose values depend on changes in the field strength. The requirement to know the limiting hysteresis loop and a few minor loops is unquestionable because a comparison between calculation and measurements results should be the basic criterion for assessing the correctness and usefulness of the proposed approximation of the hysteresis changes.
Calculation times of a single hysteresis loop are comparable. On the one hand, a comparison between these times using different models is not justified, because errors between the calculation and measurement results depend, inter alia, on the assumed number of calculation steps. On the other hand, it is more important how long the calculation times are when a hysteresis model is included in a large-dimensional system of matrix algebraic equations describing the magnetic field distribution. Often the computation time is of secondary importance as it is more important to obtain the convergence and accuracy of the solution.
Taking into account the magnetic hysteresis in the equations of the magnetic field distribution is a much more difficult issue than introducing the nonlinear but unambiguous magnetization characteristics, especially since the area of the magnetic field is divided into at least several dozen elementary segments. Very often, the finite element method is used, which is based on the differential form of the Maxwell equations. In this case, the hysteresis model is taken into account after each stable solution of the equation system of the field distribution. The equivalent reluctance network method is beneficial to introduce dependencies approximating hysteresis changes into equations of the magnetic field distribution [46,51]. This method is based on the Maxwell equations in integral form, and consequently, the components of the field strength and flux density, and not their derivatives, are unknowns in these equations. It allows easy replacement of the flux density components with the appropriate field strength components or vice versa.
If the distribution of a one-dimensional magnetic field is considered, e.g., in the crosssection of a transformer sheet, then the flux density B n in an n-th elementary segment can be expressed as follows: B n (H) = t n B nr (H n ) + (t n − 1)B nd (H n ). (18) where t n = 1 when H increases, B nr is the relationship in form of Equation (1), showing changes in the flux density for increasing H, and B nd is the relationship in form of Equation (2), when H decreases in the n-th segment.
The flux densities off all elementary segments can be written in the following form: B(H) = T. · B r (H) + (T − 1). · B d (H). (19) where T is the column vector of elements t n , H is the column vector of field strength components Hn, and symbol * denotes the multiplication of two column vectors, and the vectors B r (H), B d (H) have the form: Functions of the 'exp' type mean that each element of the column vectors k Hr (H − H 0 ) and k Hd (H 0 − H) is an argument of an exponential function.
After each stable solution of the equation system of the magnetic field distribution, it is checked how the field strength changes in individual segments relating to a given electrical sheet. If a certain component of the field strength stops increasing and starts to decrease, then the corresponding element t n of the column vector T takes the value zero. The proposed method of approximation of hysteresis changes was used in calculations of the eddy current distribution, taking into account the magnetic hysteresis in the cross-section of the transformer core [52].
It should be noted that to calculate the field distribution during rotational magnetization, vector hysteresis models should be used. An example of a hysteresis vector model is described, inter alia, in [53]. In this approach, the surface of a considered electrical steel sheet is divided into an assumed number of specified directions, and in each direction, a certain hysteresis loop is assigned. In this case, the hysteresis assigned to each direction differs from the hysteresis of the whole sample of the tested sheet.

Conclusions
The aim of the article is not to formulate a new hysteresis model that would be associated with the magnetization process on a micro scale but to propose a description of hysteresis changes that can be relatively easily introduced into the equations of the magnetic field distribution. Therefore, it seems that a comprehensive comparison of the proposed approximations with the existing hysteresis models based on the physics of this phenomenon is not substantively justified.
The proposed method, which enables the approximation of hysteresis changes in electrical steel sheets, has several advantages in relation to other hysteresis models. The changes in flux density, as a function of the relationship and inverse relationship of the field strength, are determined using simple expressions. This makes it possible to relatively easily take into account the hysteresis phenomenon in calculations of the magnetic field distribution in electrical steel sheets. Note that the described method enables satisfactory results in modelling the hysteresis changes occurring in transformer sheets for directions other than the rolling direction to be obtained, unlike existing hysteresis models.
To apply the proposed method, the limiting hysteresis loop of a particular electrical sheet should be known. Additionally, some minor loops must be measured to correctly select the attenuation coefficients of 'transient' components. In some scenarios, these coefficients can have constant values; this shortens the calculation times. It is worth emphasizing that if it is possible to assume a Bur = a Bbd = 0 in the approximation B = f (H), then only the value of the exponent p B must be determined, while for the inverse approximation H = f (B), the exponent p H should be dependent on the assumed maximum value of the flux density.
The performed numerical calculations demonstrate that the proposed method enables the approximation of the hysteresis changes for different types of electrical steel sheets with sufficient accuracy for engineering research; these sheets may have both narrow and relatively wide hysteresis loops.