Finite-Element Modeling and Characterization of Iron Losses in 12 mm Thick Steel Laminations Including the Effect of Cutting

Iron losses in laser-cut toroidal samples of 12 mm thick steel laminations used in large synchronous motors are studied. Eddy currents in the lamination cross-section are solved with the 2-D finite element method while applying a constitutive law based on the Jiles-Atherton hysteresis model. The effect of cutting on the material properties is included by a continuous local material model approach, which enables to express the material properties as a function of distance from the cutting edge. The accuracy of the model is validated by comparing the simulations and experimental measurements of five toroidal samples assembled from concentric rings with different widths. Highly accurate results are obtained in terms of both the matching of B-H loops and the total loss values with an average relative error less than 2.9%. The results show that the hysteresis loss under quasi-static excitation increases up to 20.4% due to the effect of cutting. It is observed that the eddy-current loss becomes dominant over the hysteresis loss even at 5 Hz, and this eddy-current loss decreases up to 72.5% as the number of concentric rings increases. The presented model and the results accurately show how iron losses occur in thick materials and how they are affected by the cutting process.


I. INTRODUCTION
Thick steel laminations are widely used in the rotor pole shoes of large-diameter synchronous machines due to manufacturing costs.As a result of the large lamination thickness, the skin effect becomes remarkable even at low frequencies (i.e.5-10 Hz) causing significant eddy-current loss.Usually, these materials are cut to desired shapes by different techniques for appropriate usage in the machine parts.The cutting process causes degradation in the material properties [1]- [4], which also affects the distribution of the eddy currents in the material.Thus, for accurate prediction of the iron losses in such materials, a method is needed to model the eddy currents and hysteresis when considering the effect of cutting on the magnetization and losses.
The associate editor coordinating the review of this manuscript and approving it for publication was Feiqi Deng .
In the literature, several studies have been performed to model the iron losses in magnetic materials [5]- [14].In [5], [6], the eddy currents were modeled by considering a low-frequency approach without including the skin effect.In [7]- [9], the skin effect was included by modeling the eddy currents in the lamination depth with a 1-D finite element (FE) model and coupling it to the 2-D field solution.Moreover, in [10], [11], homogenization approaches were used to account for the skin effect.With these homogenization approaches, the flux-density distribution along the lamination depth was expressed analytically by considering single-valued (SV) material properties.In [12], a similar approach was used to model the eddy currents by including the hysteretic material properties.Although these models account for the skin effect, the thickness of the sheets was negligible compared to the lamination width and the edge effects, i.e., the return path of the eddy currents at the lamination edges, were ignored.In [7], [13], [14], the edge effects were included by modeling eddy currents in 2-D in the lamination depth for comparison purposes with the 1-D eddy current approximation.Nevertheless, in these studies, the thickness of the laminations was limited to 0.5 mm, and the thick materials were not studied.
Furthermore, in order to include the effect of cutting on the magnetization and iron losses of the materials, several cutting models have been studied by researchers for the thin materials [14]- [21].In [15]- [17], the material was considered to have a degraded and a non-degraded zone, and the depth and the magnetic properties of the zones were identified based on the measurement results.On the other hand, continuous material models considering the magnetic properties of local regions as a function of distance from the cutting edge were proposed in [18]- [21].Often, the iron losses were calculated based on the analytical formulations [15]- [21].These analytical formulations are computationally fast and give reasonable results for the thin materials at low frequencies, where the skin effect is negligible.However, the analytical equations become inefficient for the thick materials or at high frequencies, due to the significance of the skin effect.Although in [14], the losses were obtained from the solutions of the field quantities rather than an analytical equation, this study was also limited to thin (i.e.0.5 mm) electrical sheets.
Based on the literature study, edge effect of the eddy currents and cutting damage in thick laminations have not been properly addressed before.In this paper, we present a 2-D axisymmetric FE model with the inclusion of a cutting model to analyze and characterize the eddy-current and hysteresis losses in the toroidal samples having thick laminations based on the experimental measurements.The eddy currents are accurately modeled in the 2-D cross-section of the samples, accounting for the edge effect.Also, the hysteresis loss is modeled based on the Jiles-Atherton (JA) hysteresis model.Inspired by the continuous local material model approach presented in [20], damaged and undamaged B-H curves are identified to define the local B-H curves as a function of the distance from the cutting edge by using a degradation profile.The presented model is simulated for several cases, and highly accurate results are obtained compared to the measurement results.

II. MEASUREMENTS A. DESCRIPTION OF SAMPLES
Measurements were performed on toroidal samples cut from typical 12 mm S275JR grade structural steel laminations used in large synchronous motors.In order to vary the amount of cutting surfaces to investigate the effect of cutting, 5 sample groups (A-E) with the same external dimensions were formed from different numbers of laser-cut concentric rings.The concentric rings were insulated from each other to avoid galvanic contacts.The geometries of the toroidal samples from the top view and cross-section are shown in Fig. 1.Specifications of these toroidal samples are given in Table 1 and Table 2.

B. MAGNETIC MEASUREMENTS
Magnetic measurements were performed between 0.25-1.5 T flux density amplitudes under the quasi-static case and sinusoidal excitations at 5 Hz and 10 Hz frequencies at ambient room temperature.The schematic of the measurement system is given in Fig. 2. The measurement system consists of a toroidal sample, an Elgar SW5250A power supply, a NI USB-6251 based data acquisition device (DAQ) connected to a PC with a MATLAB-based waveform control, and a shunt resistor for the current measurement.In order to force the flux density to be sinusoidal, the power supply is controlled by the PC through a control signal u ctrl .The average magnetic flux density B av in the toroidal sample is calculated by where u 2 is the induced back-electromotive force in the secondary winding, N 2 is the number of secondary turns, and A is the cross-sectional area of the flux path.The magnetomotive force F created by the primary winding corresponds to where N 1 is the number of primary turns and i 1 is the measured current in the primary winding.This magnetomotive force gives rise to magnetic field strength H s on the surface of the core such that which varies along the radial coordinate r.In this article, for the purpose of demonstration, we define the surface magnetic field H s,av at the mean radius as  where l av is the length of the flux path at the mean radius (see Table 1).The term B-H loop refers to dynamic relationship H s,av (B av ).The core loss density per unit mass (W/kg) during one supply period T is obtained by where ρ is the mass density.

C. ELECTRICAL CONDUCTIVITY MEASUREMENT
The electrical conductivity of the samples was measured at ambient room temperature by using the sample consisting of 1 ring (sample A).The sample was cut into 2 half pieces, and one half piece was used for the measurements.3 holes were made on both sides of the half piece and screws were placed inside the holes.The resulted geometry from the top view is given in Fig. 3.By using a DC power supply, DC currents were imposed to the sample across the arc through the screws 1-1', 2-2', and 3-3' separately, and the voltage drops were measured.After calculating the resistance R for each case from the measured voltage drop and the imposed current, the conductivity σ for each case was calculated based on the analytical solution of where d, r out , and r in are the lamination thickness, outer radius, and inner radius, respectively.The conductivity of the material was approximated as the average of the calculated conductivity values of all cases, and obtained as 5.6 MS/m.

III. MODELING A. AXISYMMETRIC FE MODEL
An axisymmetric 2-D FE model presented in [13] is used and developed further with the implementation of a JA hysteresis model and the inclusion of a cutting model.In Fig. 4, the geometry of the problem is given in rφz-coordinate system with unit vectors r, φ, and ẑ.The bold quantities represent the vectors, while the non-bold quantities represent the scalars.
The numerical analysis is based on the T formulation.The electric vector potential T and the reduced magnetic scalar potential are considered as to formulate magnetic field strength H = H φ such that Here, H s (r, t) = F(t)/2πr corresponds to magnetic field strength value at the surface of the geometry when T is set to 0 on the surface by a homogeneous Dirichlet condition.
In the axisymmetric case, Ampere's law (11) and Faraday's law (12) can be expressed as where E = E r (r, z, t)r + E z (r, z, t)ẑ, J = J r (r, z, t)r + J z (r, z, t)ẑ, and B = B(r, z, t) φ are electric field strength, electric current density, and magnetic flux density vectors, respectively.Combining ( 11)-( 12) yields F(t) from the measurements is used as the source to the problem and ( 13) is discretized by using the Galerkin FE method with test functions T such that The FE mesh consists of 800 quadratic triangular elements with 1701 nodes.The constitutive material model affects the simulated eddy-current distribution and losses.In this paper, both hysteretic and SV B(H ) relationships are studied for different analyses.Eddy-current loss coupled with hysteresis loss is simulated by using a hysteretic relationship, and uncoupled eddy-current loss is simulated by using a SV relationship.Throughout the paper, these cases are referred as coupled case and uncoupled case, respectively.
The hysteretic relationship is obtained from the scalar JA model [22].The model can be summarized by the following equations: where H and H eff are considered as applied and effective magnetic field strengths.M , M an , and M irr represent total, anhysteretic, and irreversible magnetizations, respectively.p = [α c k M s a] T is a vector of model parameters, which are obtained by fitting.Based on these parameters, the hysteretic relationship is defined as B hy (p, H ). The SV relationship B sv (p, H ) is obtained from the modeled B hy (p, H ) by using the JA model to simulate the minor loops at different flux density amplitudes and then tracing the peak points of these minor loops.
After obtaining the solution, the instantaneous rate-ofchange of the magnetic field energy p hy (t) and eddy-current loss p cl (t) averaged over the volume V are obtained by In order to obtain the time-averaged hysteresis p hy and eddy-current losses p cl , p hy (t) and p cl (t) are averaged over one period T of a closed cycle such that

B. MODELING THE EFFECT OF CUTTING
The degradation caused by the laser cutting process is modeled by following the continuous local material model approach presented in [20].The model expresses the magnetic flux density B as a function of the magnetic field strength H and distance from the cutting edge x by using a magnetization curve for the undamaged material B un (H ), a magnetization curve for the damaged material B dam (H ), and a degradation profile η(x) by In this work, a quadratic degradation profile given by is used, where δ is the degradation depth, beyond which the degradation disappears.
Based on the formulations in ( 24)-( 25), the undamaged curve B un (H ), the damaged curve B dam (H ), and δ need to be identified for the hysteretic and SV cases.These two curves cannot be measured directly.Instead

A. IDENTIFICATION OF PARAMETERS
The identification procedure of the JA model parameters p un and p dam and the degradation depth δ is discussed in Sec.III-B.By following the steps for the identification procedure, initially, the hysteretic undamaged curve B un,hy (H ) is identified by least-squares fitting by comparing the simulated and measured B-H loops of sample A at 0.5 T, 1 T, and 1.5 T amplitudes for the quasi-static case, 5 Hz, and 10 Hz frequencies simultaneously.For the simulations of the quasistatic case, the frequency is used as 0.001 Hz.After identifying the hysteretic undamaged curve B un,hy (H ), the hysteretic damaged curve B dam,hy (H ) and the degradation depth δ are identified similarly by least-squares fitting by comparing the simulated and measured B-H loops at 0.5 T, 1 T, and 1.5 T amplitudes for the quasi-static case of all samples simultaneously.While fitting the JA parameters p un and p dam , M s is kept same for both cases as it is a material property and not affected by cutting process.As a result of fitting, the degradation depth δ is obtained as 4.1 mm.The fitted JA parameters for the identified hysteretic B-H curves are given in Table 3.The identified hysteretic B-H curves based on the fitted JA model parameters and the degradation profile are given in Fig. 5(a).The identified SV undamaged and damaged curves based on the simulated minor loops of the identified hysteretic curves are given in Fig. 5(b) and Fig. 5(c).

B. SIMULATION RESULTS
By using the identified hysteretic B-H curves and the degradation depth, the total losses for the coupled case are simulated for all toroidal samples at different flux density amplitudes and frequencies.The simulated and measured B-H loops of samples A, C, and E at 0.5 T, 1 T, and 1.5 T amplitudes are presented in Fig. 6 for the quasi-static case, 5 Hz, and 10 Hz frequencies.It should be noted that in the quasi-static case, the eddy currents do not exist and all the losses are due to the hysteresis loss.Fig. 6 shows that the simulated and measured B-H loops match successfully.Although the degree of freedom is high while fitting the model parameters by considering different measurements of different samples simultaneously, the simulated results show that identified model parameters yield consistent results for all cases, which shows the accuracy of the identification procedure.The results for the quasi-static case in Fig. 6 show that the average permeability of the material is significantly affected by the cutting, which also affects the hysteresis loss of the material.When the frequency is 5 Hz, the area of the B-H loops increases remarkably compared to the quasi-static case as a result of the significance of the eddy-current loss.This significance increases further as the frequency increases to 10 Hz.
Similar to the simulations of the coupled case, the simulations are performed for the uncoupled case by using the identified degradation depth and the SV undamaged and damaged B-H curves obtained from the identified hysteretic B-H curves.The losses obtained from measurements and simulations for all cases are presented and compared in the next section in detail.

C. COMPARISON OF THE LOSSES
The simulated and measured losses for all cases are given in Fig. 7.The simulations of the coupled case yield the simulated total losses, which are segregated into the hysteresis and coupled eddy-current losses.The simulations of the uncoupled case yield the uncoupled eddy-current loss.
The results for the quasi-static case in Fig. 7 show that the simulations match with measurements with the average relative error less than 5.4% for all simulations.These losses correspond to the hysteresis loss as the eddy currents do not exist in the quasi-static case.It is seen that the hysteresis loss increases as the number of concentric rings increases due to the increasing deformation introduced by the cutting.For instance, at 1.5 T, the hysteresis loss of sample E is 20.4% larger than the hysteresis loss of sample A.
Similarly, the results for 5 Hz and 10 Hz in Fig. 7 show that the simulated total losses match accurately with the measured total losses for all samples with the average relative error less than 2.9% for all simulations.The comparison of the segregated loss components indicates that the major part of these losses is composed of eddy-current loss even though the studied frequencies are low.This is due to the fact that the studied material is considerably thick, and the skin effect is remarkable, which results in a large amount of eddy-current loss.
The significance of the eddy-current loss in proportion to the simulated total losses increases as the frequency increases.For instance, the coupled eddy-current loss for sample A is 83.7% of the simulated total loss at 5 Hz, 1.5 T, while this proportion increases to 91.1% at 10 Hz, 1.5 T. On the other hand, under the same frequency level, the significance of the eddy-current loss decreases as the number of concentric rings increases.As the number of concentric rings increases, the electromotive force which gives rise to the eddy  currents in each ring decreases.At the same time, the overall resistance experienced by the eddy-current loop in the ring cross-section remains roughly constant since the resistance along the radial direction decreases, while the resistance along the thickness increases.Therefore, the reduced electromotive force with the approximately constant overall resistance results in a lower current density, and the eddy-current loss reduces.The distribution of current density for the samples A and E at 5 Hz, 1.5 T is given in Fig. 8 for illustration.A 1-D eddy-current loss model would assume the sheet to be infinitely long in the radial direction, making it impossible to account for the above-mentioned effects.
Moreover, as the number of concentric rings increases, a similar decreasing trend is observed in the simulated total losses.This is due to the decrease in eddy-current loss as the hysteresis loss does not change significantly compared to the eddy-current loss.For instance, the comparison of the losses of sample A with the sample E at 5 Hz, 1.5 T shows that although the hysteresis loss increases from 2 W/kg to 2.3 W/kg, the coupled eddy-current loss decreases from 10.2 W/kg to 3.3 W/kg, which decreases the total losses from 12.2 W/kg to 5.6 W/kg.Also, the comparison between the coupled and uncoupled eddy-current losses in Fig. 7 shows that, there are differences between their values and behaviors at different flux density amplitudes.The differences imply that the eddy current distribution in the materials is slightly affected by the hysteretic characteristic of the materials.These results are in line with the findings of [23], where the dependence of the eddy currents on hysteresis was studied.

V. CONCLUSION
In this paper, we proposed a combined experimental and numerical procedure to characterize material properties and 115716 VOLUME 9, 2021 iron losses of very thick steel laminations.The method removes the need of field homogeneity, which would be impossible to achieve in thick samples.For accurate prediction of the iron losses, the eddy currents should be modeled in 2-D to account for edge effects correctly.The results show that the interdependence between the eddy current and hysteresis losses are affected by the level of magnetization, but the effect of this interdependency on the total iron losses still need to be investigated.

FIGURE 1 .
FIGURE 1. Geometry of the toroidal samples from the top view and cross-section.

FIGURE 2 .
FIGURE 2. Schematic of the measurement setup.

FIGURE 3 .
FIGURE 3. Top view of the constructed geometry for the electrical conductivity measurement.There are 3 holes on the both sides and screws are placed inside the holes.

FIGURE 4 .
FIGURE 4. Geometry of the problem for the axisymmetric FE model.

p
dam that are identified in an iterative manner by comparing simulation results against measurements.The identification procedure is implemented by the following steps: 1) Assume that sample A is undamaged.Run a least-squares algorithm to fit the JA parameters p un for the hysteretic undamaged curve B un,hy (H ) = B hy (p un , H ) such that the FE-simulated B-H loops match with the ones measured from sample A at different frequencies and amplitudes.2) Using the hysteretic undamaged curve B un,hy (H ) found at step 1, run a least-squares algorithm to fit the JA parameters p dam for the hysteretic damaged curve B dam,hy (H ) = B hy (p dam , H ) as well as the degradation depth δ such that the FE-simulated B-H loops match with the ones measured from samples B, C, D, and E at different frequencies and amplitudes.3) Obtain SV undamaged B un,sv (H ) = B sv (p un , H ) and damaged B dam,sv (H ) = B sv (p dam , H ) curves from the identified hysteretic B un,hy (H ) and B dam,hy (H ) curves.

FIGURE 5 .
FIGURE 5. Identified hysteretic B-H curves and degradation profile (a), identified SV undamaged curve (b), and identified SV damaged curve (c).The undamaged and damaged curves meet at full saturation.

FIGURE 6 .
FIGURE 6. Simulated and measured B-H loops of samples A, C, and E at 0.5 T, 1 T, and 1.5 T amplitudes.Note the difference in the scales of x-axis in the graphs.

FIGURE 7 .
FIGURE 7. Simulated and measured losses of samples for all cases.

FIGURE 8 .
FIGURE 8. Rms current density distribution and isolines of the rms field strength in (a) sample A and (b) sample E at 5 Hz, 1.5 T.

TABLE 1 .
Specifications of the toroidal samples I.

TABLE 2 .
Specifications of the toroidal samples II.
, both curves are described by the JA model with different parameters p un and

TABLE 3 .
Fitted JA parameters of the identified hysteretic B-H curves.