Finite element level validation of an anisotropic hysteresis model for non-oriented electrical steel sheets

This paper presents the finite element level validation of the anisotropic Jiles–Atherton hysteresis model. Numerical analysis of a round rotational single sheet tester is performed using the 2D finite element method. Anisotropic extension of the Jiles–Atherton hysteresis model is coupled with the 2D finite element method. The finite element simulations are performed for the cases when the magnetic field alternates and rotates in the lamination plane. The simulated results for alternating and rotational flux density excitations agree with the measured data. The measured data used in this paper corresponds to the M400-50A nonoriented silicon steel.


Introduction
Nonoriented (NO) silicon steels are suitable for electromagnetic applications where the core experience directionally varying magnetic fields, for example, the stator yoke of a medium-sized induction motor [1].The magnetic field can alternate and rotate in the stator core, depending on the location.A phenomenological hysteresis model should predict the magnetic fields for different types of excitations: alternating, rotating, and elliptical.The model should also describe magnetic responses under external stress and thermal loadings [2][3][4].However, a physics-based hysteresis model of such caliber is unavailable [5].
The NO silicon steel possesses a significant level of magnetic anisotropy [6][7][8][9].Several models of magnetic anisotropy have been proposed for NO silicon steel.In [10], a method based on magnetic energy is proposed to account for magnetic anisotropy observed in the NO silicon steel.The model expresses energy with a Gumbel-type distribution, whose parameters are extracted utilizing the anhysteretic data.The model is suitable for the finite element method (FEM).However, the results they predict disagree with the measurement data.Instead of using the magnetic energy, [11] employs bicubic splines to represent the anhysteretic characteristics identified from the rotational measurements.The method based on the surface bicubic spline is easier to integrate with the FEM and it is more accurate than the energy-density method.Contrarily, [12] expresses magnetic energy using the Fourier series.The coefficients of the Fourier series-based model are identified using anhysteretic data in several measurement directions.[13] shows that the Fourier series-based anisotropic model is well adapted to the FEM by simulating the magnetic field in an interior permanent magnet machine composed of NO silicon steel.
Alternatively, [14] has modified Mayergoyz's approach to model anisotropy observed in a NO silicon steel sheet.The modification adds the Everett data from the alternating measurements in several directions.The anisotropic generalization of Mayergoyz's model controls the amplitude of the input projections, introducing anisotropy in the predicted results.[15,16] extend vector version of the Jiles-Atherton (JA) model to include magnetic anisotropy.Their extension considers measurement data in only two directions-the rolling and transverse.The anisotropic extension of the JA model estimates the field strength loci based on the parameters from rolling direction (RD) and transverse direction (TD) in the intermediate directions.Likewise, [17,18] implement anisotropy in the energy-based (EB) hysteresis model to predict variations in alternating and rotational field strength in a NO silicon steel sheet.The anisotropic EB model utilizes parameters identified from the alternating measurements in the RD and TD.Moreover, the pinning field probability density is represented by a univariate Gaussian-type distribution, whose parameters are estimated by fitting the model to the measurement data [17].
Likewise, [19] presents an anisotropic play-type model similar to [17]'s.The model considers parameters from the alternating measurements in all principal axes: the RD, TD, and the silicon steel lamination's thickness.The anisotropic play-type model shows a reasonably good fit for the magnetic field alternating in the RD, TD, https://doi.org/and 45 • directions.Moreover, [20] proposes an anisotropic play-type model.The model introduces anisotropy in the isotropic play-type model with the parameters identified from several unidirectional alternating measurements.The anisotropic play-type model yields a reasonably good fit with the measured data for NO silicon steel.
A physics-based multi-scale model (MSM) can be employed to predict the anisotropic behavior in a NO silicon steel lamination [21].However, the MSM uses texture data, requiring complex measurements.Additionally, the MSM is a large minimization problem and computationally heavy.For the FEA, a simplified MSM is considered in [22].The anhysteretic magnetization identified from the simplified MSM is fed to the hysteresis models-such as JA, EB, and Hauser [23,24].In [22], the simplified MSM-JA model is coupled with the FEM to study the effect of stress-induced anisotropy in a switched reluctance motor.Moreover, [25] couples the MSM with the EB model to account for the anisotropy observed in grain-oriented silicon steel.Nevertheless, the physics-based MSM-EB or MSM-JA model could be a good choice for predicting magnetic behavior in the NO silicon steel [22,25].
Characterization of NO silicon steel is an essential part of the electrical machine design process [26].An Epstein frame is a standardized device for measuring alternating BH characteristics.However, the Epstein frame is limited to alternating fields.A rotational single sheet tester is employed to measure rotational fields because it allows measurements of vector components using the field sensing methods [27].A round rotational single sheet tester (RRSST) is a popular choice as it ensures better rotation of the flux density, utilizing a larger measuring region.Additionally, the RRSST allows measurement of alternating and rotational magnetic flux density-magnitudes as high as 2 T [28,29].
This paper presents the finite element analysis (FEA) of an RRSST sample-a NO silicon steel disk of grade M400-50 A. The FEA considers anisotropic and isotropic Jiles-Atherton (JA) hysteresis models [30].The following four cases are considered for the FEA: flux density alternating in the RD, alternating in the 45 • direction (ALT45), alternating in the TD, and rotating in the counter clock-wise direction (ROT).The anisotropic JA model (AJAM) utilizes the anhysteretic magnetization from alternating BH measurements in seven directions: 0 • , 15 • , 30 • , … , 90 • [30].Moreover, the JA model's parameters vary as a function of the magnitude and polar direction of the flux density [31].The unidirectional alternating BH measurements are used to estimate the JA model's parameters.In contrast, rotational measurements are used for validation (see Table 1).

Anisotropic JA model
This paper considers AJAM explained in [30]-an improved extension of Bergqvist's vector JA model [32].The model equations are re-summarized in the following:

Differential susceptibility
Bergqvist's vector JA model defines differential susceptibility [15]: where  is the identity matrix,  and  are the differential anhysteretic and irreversible susceptibilities,  ≥ 0 is a parameter that represents interdomain coupling, and  ∈ [0, 1] is another parameter that quantifies the reversible processes associated with the bowing of domain walls ( = 1 corresponds to a completely reversible process) [33,34].The differential susceptibility in (1) is expressed using an auxiliary vector as follows: where  ≥ 0 is a parameter associated with the pinning of domain walls,  an represents the anhysteretic magnetization, and  ef f =  +  is the effective field strength.The condition in (2) reflects the assumption that irreversible changes occur only when the field  ef f is incremented in the direction of (  an −  ) .

Differential reluctivity
The FEM based on magnetic vector potential formulation uses differential reluctivity [35]: where  0 is the permeability of the free space,   represents the differential permeability, and   is the differential susceptibility given by (1).

Anhysteretic magnetization
The anhysteretic magnetization for the isotropic case is usually expressed as [35] where  ef f = ‖ ef f ‖ is the norm of effective field, and  an =  (  ef f ) describes the anhysteretic curve identified from the unidirectional alternating BH loop.The differential anhysteretic susceptibility is expressed as In this paper,  an defined in ( 4) is identified by averaging the unidirectional alternating BH characteristics measured in seven directions: Because the measurement is performed at an excitation frequency of 50 Hz, the hysteretic field strength is extracted by subtracting the eddy-current loss field from the measured field strength  meas : where  and  represent the electrical conductivity and thickness of the NO silicon steel sample.The derivative on the right-hand side of ( 6) is evaluated using the central difference method.After that, the phase shift between  and  is neglected, meaning  and  are projected in the reference direction, which is the polar direction of vector .This process is performed for all the measurement directions.The field strength is averaged because the measurement is B-controlled [31]: where   = ( − 1)15 • is the polar direction of ,  = 1, 2, 3, … ,  represents the data point index, and  is the total number of data points in one excitation cycle.Thus,  an is obtained by averaging the ascending and descending branches of the major ( avg )-loop [36].
The anhysteretic magnetization follows a more generic expression for the anisotropic case: where  is the polar angle of  ef f .The components of the anhysteretic magnetization (8) could be modeled using the multiscale approach [21].However, the multiscale approach is a large minimization problem and computationally heavy.Analytical representation of the anhysteretic characteristic is possible with the help of transcendental functions [25].Nevertheless, to ease the computational burden, this paper applies bicubic splines [30].Explicitly, in each rectangle k, where  k x and  k y are 4 × 4 coefficient matrices whose elements are identified from the anhysteretic M(H) characteristics (see [30,Fig. 2]).The source code available in [37] is used to obtain the elements of  k x and  k y .The elements of  are evaluated from ( 9) and ( 10) using the following polar-to-Cartesian transformation:

Directional variation of JA parameters
The following equations describe the directional variations of the JA model's parameters [31]: where  JA = {, , },  x = { x ,  x ,  x },  y = { y ,  y ,  y }, and  and  represent the norm and polar direction of .The parameters are identified from the unidirectional alternating BH measurements in only two directions-the RD and TD.The parameters in the right hand side of ( 12) are represented with a piece-wise linear polynomial along the magnitude of the flux density (see [30,Fig. 5]).

Finite element modeling
This section considers a magnetostatic field problem in a domain  ∈ R 2 .The domain is a collection of linear  L and nonlinear  N subdomains:  =  L +  N .The spatial variation of the magnetic field strength in  is described by the following partial differential equation: where  s =  (, ) z is the source current density-usually imposed on the subdomain  L .A hysteretic material law  = ℋ () links  and  in  N , whereas, in  L , a linear relationship  =  0 , where  0 = 1∕ 0 , exists between them.
For any continuous vector potential  = (, ) z , the flux density  = ∇ × , and  satisfies the Coulomb gauge, ∇ ⋅  = 0.The Galerkin method is employed to solve (13).The equation is multiplied by a weight function:  =  z , and integrated over the domain , resulting in [38] where (, )  and ⟨, ⟩  denote the integration ∫  ( ⋅ )d and ∮  ( ⋅ )d , and  is the outward oriented unit vector normal to the boundary  of .For simplicity, the contour integral on the right hand side of ( 14) is omitted.The boundary conditions will be discussed in a separate section.

Solution of the nonlinear field equation by fixed point method
The following constitutive relation is used to solve ( 14) by the fixed-point method [39,40] where  FP is a reluctivity like quantity, and  is a magnetization like quantity.Thus, using (3), Moreover, using ( 15) in ( 14) results in the linearized system: In the FEM, the approximation of the vector potential where   and   are the nodal value and shape function associated with node  of the finite element (FE) mesh, and  represents the total number of nodes.Commonly, the shape functions connected to the free nodes are used as the weight functions.This paper adopts the triangular elements to discretize , utilizing linear shape functions so that the weight functions are linear basis functions [38].
Eq. ( 17) is expressed in the following form: where  = [ where  =  N +  L is the system matrix, which is usually sparse and symmetric, and   +1 =  +1 −  +1 is the load vector, subscripts +1 and  denote the current and previous time instant, and superscripts  + 1 and  denote the current and previous iterate.This paper uses the sparse direct linear solver to solve the system of linear Eqs.(20) [41].The nonlinear iterations in (20) are stopped if the  1 norm of the change in solution ensures where  abs and  rel represent the absolute and relative tolerances, and  min and  max denote the minimum and maximum nonlinear iterations.

Boundary condition
The following boundary condition is applied on the outer boundary of the RRSST sample to set the magnetic flux density to alternate in any arbitrary direction : where  b is the -component of the vector potential on the outer boundary  , Âb is the peak amplitude of the vector potential,  is the angular coordinate of the particular point on the boundary,  = 2 is the angular frequency, and  is the frequency in cycles per second.Moreover, the flux density rotates on the sample if the vector potential satisfies the following condition [42]: The flux density's amplitude B is controlled by controlling the value of Âb .Thus, using Âb =  B, where  is the radius of the circular steel sheet, in ( 22) and ( 23), respectively.Furthermore, considering the Dirichlet boundary condition, ( 20) is expressed as follows: where the additional subscripts h and b denote the free and fixed boundary nodes.

FE mesh of the RRSST
The layout of the RRSST is detailed in [30,Fig.4].Fig. 1 is the 2D FE mesh of the circular steel sample.The mesh is generated using the freely available meshing software Gmsh [43].The mesh consists of two material domains: air and silicon steel.The elements with yellow edges represent the air domain, whereas those with red and blue edges represent silicon steel.The elements with blue edges distinguish the sensor region from the remaining core.The mesh consists of 29,682 elements, of which 144 are air elements, and the remaining 29,538 elements represent silicon steel core.There are 14,898 nodes, of which 112 belong to the outer boundary, and the remaining 14,786 are internal free nodes.

Estimation of 𝑩 and 𝑯
The magnetic flux density in the sensor region (see Fig. 1) is obtained by using the following expression: where denotes the average value of the -component of the vector potential,   is the cross-sectional area of the  th hole, and  w is center-to-center distance between holes aligned in the RD or TD.Note that the holes with index  = 1 and 2 are aligned in the TD, and holes with index  = 3 and 4 are aligned in the RD.Likewise, the components of the magnetic field strength sensed by the H-coils, which is placed on the surface of the sample, are considered to be the weighted average: where  is the total surface area of the sensor region excluding the holes.Moreover, after estimating  and  from ( 25) and ( 26), the hysteresis loss density is evaluated using the following integral [44]: The integral in the right hand side of ( 27) is numerically evaluated using the trapezoidal rule.

Evaluation of JA parameters at Gauss integration points
The surface integral in the 2D FEM is obtained using the Gauss-Legendre quadrature rule [38].The surface integration is converted to a weighted sum according to the quadrature rule.The parameters of the JA model are evaluated from the flux densities at the two previous steps at the integration point (IP): where  JA = {, , }, and subscript  denotes the time step.This paper considers only one IP-the centroid of the reference triangle.

The R-squared measure of goodness of fit
The R-squared ( 2 ) is used as a quantitative measure of goodness of fit between the predicted and measured magnetic data.The following equation defines the  2 : where  meas and  simu denote the measured and predicted values, Ȳ is the mean value, and  denotes the total number of data points.For the RRSST, the goodness of fit is evaluated by setting  = { x ,  x ,  y ,  y } in (29).For an exact match between the predicted and measured results,  2 = 1.Thus, a positive value of  2 close to 1 could be considered a good fit.

Results and discussion
This section presents the results of the 2D FEA of the RRSST sample.Eqs. ( 22) and ( 23) are considered to establish either unidirectionally alternating or rotating magnetic flux density in the circular steel sample.The magnitude of  is controlled by setting B = {0.5, 1, 1.5}  and  = 1 Hz in the equations describing the essential boundary conditions.The FE simulations are performed for the first five periods of the supply frequency.A single period of the supply frequency is discretized into 400 steps.The number of minimum allowed nonlinear iteration is set to be  min = 2.The values of tolerances are fixed as  abs = 10 −5 and  rel = 10 −7 .Equations ( 25) and ( 26) are used to estimate  and  in the sensor region of the RRSST sample.The last or fifth period of the magnetic field solution is considered to obtain the hysteresis loss density.The results of the FE simulations are in Figs.2-10.The following subsections discuss the results.

𝑩 Alternating in the RD
Fig. 2 shows the result when the magnetic flux density alternates parallel to the RD.The isotropic JA model (IJAM) considers BH characteristics averaged over seven directions, meaning the prediction would be closer to the magnetic characteristics in the 45 • direction.Because the field strength is higher for the isotropic case, the losses would also be higher for the isotropic case than anisotropic.Fig. 3 compares the measured and simulated alternating BH characteristics in the RD.The results of the AJAM agree with the measured data.However, some differences are seen in the results of the IJAM.Clearly, the IJAM underestimates the BH characteristics when the flux density alternates in the RD.

𝑩 Alternating in the 45 • direction
Fig. 4 shows the result when the magnetic flux density alternates in the 45 • direction.The differences between the results are difficult to see from the field distribution.The hysteresis loss distribution shows noticeable differences.The losses are higher for the AJAM than for the IJAM.Fig. 5 compares the measured and simulated BH loci.The BH characteristics predicted by the AJAM agree better with the measurement data than the IJAM.However, notable differences are seen at a high magnitude of the flux density.

𝑩 Alternating in the TD
Fig. 6 shows the result when the magnetic flux density alternates in the TD.Because the TD is a hard direction of magnetization, the magnetic field strength is higher for the AJAM than IJAM.Consequently, the losses are higher for the anisotropic case.Fig. 7 compares the measured and simulated alternating BH characteristics.The AJAM's results are in good agreement with the measured data.In contrast, the IJAM's results differ significantly from the measured data.Thus, the IJAM underestimates the BH characteristics if the flux density alternates in the TD.

𝑩 Rotating in the counter-clockwise direction
Fig. 8 shows the result when the magnetic flux density rotates in the counter clock-wise direction.The field quantities rotates at each step, so, it may not be intuitive to compare the results from two different models at a single time instant.Nevertheless, the hysteresis losses are evaluated from a full cycle of the input excitation, so, comparing the isotropic and anisotropic models is more intuitive.The distribution of the hysteresis loss density is more pronounced for the AJAM than IJAM (see Figs.  of the AJAM agree with the measured data until 1 T.However, the differences are significant at 1.5 T. Contrarily, the IJAM's results differs considerably from the measured data.

Comparison of R-squared coefficient
Tables 2 and 3 show the  2 values for the FEA results of the IJAM and AJAM.Note that the  2 values are tabulated for the components of .For the components of ,  2 = 0.9999 in all cases, so it is not tabulated.The  2 values for the AJAM are higher than those for the IJAM, meaning the predictions are better for the AJAM.The  2 values for IJAM decrease with an increase in excitation magnitude, signifying the predictions worsen at high flux density excitations.

Comparison of hysteresis losses
Fig. 10 compares the hysteresis losses for four different types of input excitations.The IJAM should yield identical magnetic characteristics for alternating input excitations, which means the losses should also be the same.Thus, according to the result presented in Fig. 10(a), input excitations along RD, ALT45, and TD yield similar losses.In contrast, for the anisotropic model, alternating losses depend on the direction of the input excitation.As shown in Fig. 10(b), losses produced by the field alternating in the RD is the lowest, followed by ALT45 and TD.If we now compare rotational losses, we see that the differences are negligible until 1 T, and at 1.5 T, the AJAM yields higher losses in the RRSST sample.

Simulation time and iterations
The FE simulations are performed using Intel ® Core™ M-5Y51 @1.10 GHz RAM 8 GB with GNU/Linux environment.An in-house program written in the C programming language is used for the simulations.The information related to the simulation time and average number of nonlinear iterations are in Tables 4 and 5.The FE simulation that considers the IJAM spends about one second per step on average, whereas, the simulation times are on average 10% higher for the AJAM.Nevertheless, the average numbers of nonlinear iterations are the same.

Fig. 1 .
Fig.1.2D FE mesh of the RRSST sample.The sensor region includes four holes that are drilled to accommodate the B-coils.The radius of the circular steel sample  = 39 mm, and the hole  hole = 0.4 mm (see[30, Fig.4]).A zoomed section on the right hand side depicts the triangular mesh inside the hole that is aligned in the RD.

Fig. 2 .
Fig. 2. Distribution of the magnetic field strength and flux density at a time instant  = 4.25 s, and time-averaged hysteresis losses in the RRSST sample.The magnetic flux density alternates in the RD.(a), (c) and (e) FEM coupled with the IJAM.(b), (d) and (f) FEM coupled with the AJAM.

Fig. 3 .
Fig. 3. Measured and FE simulated BH loops captured in the sensor region of the RRSST sample.The magnetic flux density alternates in the RD.(a) FEM coupled with the IJAM.(b) FEM coupled with the AJAM.

Fig. 4 .
Fig. 4. Distribution of the magnetic field strength and flux density at a time instant  = 4.25 s, and time-averaged hysteresis losses in the RRSST sample.The magnetic flux density alternates in the 45 • direction.(a), (c) and (e) FEM coupled with the IJAM.(b), (d) and (f) FEM coupled with the AJAM.
8(a) and 8(b)).The presence of easy and hard magnetization directions leads to different amounts of losses in different directions around the perimeter of the holes.The rest of the sample experiences uniform loss density distribution.Fig.9compares the  x ( y ),  x ( y ),  x ( x ), and  y ( y ) characteristics.The results predicted by the AJAM agree with the measured data.However, notable differences are seen between the measured and simulated data at a high magnitude level.Conversely, the IJAM's prediction differs significantly from the measured data.Figs.9(e)-9(h) compare the measured  x ( x ) and  y ( y ) characteristics.The results B.Upadhaya et al.

Fig. 5 .
Fig. 5. Measured and FE simulated  x ( x ) and  y ( y ) loci obtained from the sensor region of the RRSST sample.The magnetic flux density alternates in the 45 • direction.(a) and (c) FEM coupled with the IJAM.(b) and (d) FEM coupled with the AJAM.

Fig. 6 .
Fig. 6.Distribution of the magnetic field strength and flux density at time instant  = 4.25 s, and time-averaged hysteresis losses in the RRSST sample.The magnetic flux density alternates in the TD.(a), (c) and (e) FEM coupled with the IJAM.(b), (d) and (f) FEM coupled with the AJAM.

Fig. 10 .
Fig. 10.FE simulated hysteresis losses in the RRSST sample for the cases when the magnetic flux density alternates in the RD, 45 • , TD, and rotates in the counter clock-wise direction.(a) FEM coupled with the IJAM.(b) FEM coupled with the AJAM.

Table 1
Magnetic measurements for the parameter identification and model validation.

Table 4
FE simulation time and nonlinear iterations for the IJAM.