Iron Losses in Electromagnetic Devices : Nonlinear Adaptive MEC & Dynamic Hysteresis Model

In this paper, an original approach allowing the determination of the iron losses in the electromagnetic devices is presented. This new approach exploits the Loss Surface (LS) hysteresis model and the magnetic flux density waveforms resulting from a generalized nonlinear adaptive magnetic equivalent circuit (MEC) using a mesh-based formulation in two-dimensional (2-D) or quasi three-dimensional (3-D). The model coupling has been applied to a 18-slots/16-poles radial-flux interior permanent-magnet (PM) synchronous machine (PMSM) dedicated to automotive applications, mainly for electric/hybrid/fuel cell vehicles (EVs/HEVs/FCVs). The obtained results have been compared with those made retrospectively in the 2-D transient finite-element (FE) FluxTM. The influence of the MEC discretization on the iron loss calculation and the electromagnetic performances has been analyzed. The computation time is divided by 3/2 with an error less than 7 %.


Context
Electrical steels are highly engineered to offer a combination of the required physical and mechanical properties to successfully operate in the electrical machine while providing optimum energy efficiency.The design and optimization process of such applications require precise and reliable tools for the prediction of iron losses which is one of the pertinent properties to achieve high levels in terms of energy efficiency.
In this framework, the designers have a wide choice of models and appraoches to describe the iron losses.Among these models, we meet those based on the more general equation of Steinmetz [1].Since the Steinmetz discovery, many studies have contributed to improve the first representation of Steinmetz taking into account the unconventional waveforms.They have led to the: • Modified Steinmetz equation (MSE) [2]; • Generalized Steinmetz equation (GSE) [3]; • Improved generalized Steinmetz equation (iGSE) [4]; • Natural Steinmetz extension (NSE) [5].
A second models category brings together the formalisms based on the magnetic loss separation assumption which owes its popularity to Bertotti [6][7][8] who separates the iron losses into three contributions: hysteresis losses, classical eddy current loss and excess losses.The latter having been established first for sinusoidal waveforms, it has been improved in [9] for non-sinusoidal and rotational magnetic flux density waveforms.To those categories, are added the hysteresis loss models.Among these models, we cite: • Preisach hysteresis model in its both static [10] and dynamic [11], [12] versions; • Opera hysteresis model [13], [14] available in the commercial Opera simulation software of COBHAM society (from UK); • Viscosity-based magnetodynamic model [15]; • So-called improved vector play model [16] which recently has been successfully implemented in 2-D and 3-D transient FE solvers of ANSYS corporation (e.g., Maxwell software).
In the same way, in Grenoble, the method investigated during the last decade is the so-called Loss Surface (LS) [17][18][19].Nowadays, the iron loss prediction using the LS dynamic model of hysteresis is implemented in the FE Flux TM software of CEDRAT S.A. society (recently acquired by the american society ALTAIR) as a post-processing module for a number of common non-oriented electrical steels.All these categories of methods/models differ in terms of: • Ability to consider complex induction waveforms (e.g., including higher harmonics); • Ability to consider rotating fields; • Introduction of physical considerations; • Number of experimental data for modeling; • Computation time as well as accuracy.
The models using the Steinmetz equation together with the loss separation models are usually applicable in the literature due to their ease of use and simplicity of implementation.Furthermore, such models can only provide an indication of the order of feature of the iron losses that are not always either well defined or well quantified and often calculated without the incorporation of the additional phenomena of dissipation.Indeed, in ac regime, complex dynamic phenomena are placed on electromagnetic systems.For example, in the electrical machines, such phenomena can be induced by the geometry of the magnetic circuit, the PM magnetization pattern as well as the shapes, the magnetic leakage...These phenomena affect the material, viz., the magnetic flux density B and field H distributions in the magnetic circuit in terms of waveform and magnitude.So, it is conveniently to resort to more complex loss models, i.e., the behavioral models of hysteresis to enable a more precise prediction of iron losses, especially in the field of electrical machines, where the precision of computation is of prime importance.The loss models require, as an input, the results of the electromagnetic simulation of the device, whose purpose is to precisely determine the evolution of the B(t) signal in the magnetic circuit.To do so, the electromagnetic modeling must incorporate, depending on the investigated problem, the different encountered phenomena: the rotor motion (for the electrical machines), the coupling with the electrical/mechanical/thermal circuit equations, the conductors overheating, the power supply (i.e., the nature and the waveforms), the end-effects...The modeling techniques can be classified in various categories [20]: • Numerical methods (e.g., FE [21], finite-difference [22], boundary-element analysis [23]); • Electrical/Thermal/Magnetic equivalent circuit (EEC/TEC/MEC) [24,25]; • Schwarz-Christoffel (SC) mapping method [26]; • Maxwell-Fourier methods [27][28][29]: i) Multi-layers models, ii) Eigenvalues model, and iii) Subdomain technique.
The numerical models, even if they seem to be irreproachable in terms of precision, they are not often very suitable especially for optimization procedures because of the long computation time.Thus, the designer are often reoriented to the way of (semi-)analytical models 1 (i.e., EEC/TEC/MEC, 1 This type of model consists of N interrelated analytical equations which must be solved numerically (a.k.a.semi-numerical models).For example, in Maxwell-Fourier methods, the unknown coefficients of the series are computed by solving a (non)linear matrix system.
the SC mapping and Maxwell-Fourier methods).In comparison with the numerical methods, these models have the advantage to be explicit, accurate anf fast.Moreover, they allow us to rigorously take into account the geometry of the domains and their (non-)conductive proprieties.Taking into account the saturation effect in the SC mapping and Maxwell-Fourier methods is still a key challenge which is rarely explored in the literature.It is interesting to note that recently new scientific contributions based on Maxwell-Fourier methods have been developed considering the saturation effect [20,[30][31][32].Dubas et al. (2016) [20] realized an overview on the existing (semi-)analytical models in Maxwell-Fourier mehtods with the local/global saturation.For decades, the saturation effect can be overcame by using a nonlinear MEC.In this type of modeling, we notice that the iron losses are estimated with the aid of models such as Steinmetz [33][34][35][36][37][38], Bertotti [39][40][41][42][43] and others [44], established for the magnetic material constituting the device.Nowadays, there are no bibliographic data available in the literature dealing with the evaluation of the iron losses by using a precise hysteresis model based on the magnetic flux density waveforms resulting from a MEC.

Objectives
In this paper, and in the fast design context of the electromagnetic devices, we propose an original method for the evaluation of the iron losses in these systems.This new approach exploits the LS hysteresis model and the magnetic flux density waveforms resulting from a 2-D/quasi 3-D generalized nonlinear adaptive MEC using a mesh-based formulation.
Section 2 presents the development of the 2-D/quasi 3-D generalized nonlinear adaptive MEC allowing the computation of the magnetic state of the electromagnetic devices (e.g., in the electrical machines).The general assumptions, the automatic mesh, the matrix formulation, and the solving of the semi-analytical model have been developed in this section.The description of the LS model, the modeling of the non-oriented SiFe steel constituting the studied machine (based on experimental characterizations) as well as the implementation of the model in the FE Flux TM software have been detailed in Section 3. Section 4 describes the coupling established between the 2-D/quasi 3-D generalized nonlinear adaptive MEC and the LS model.In Section 5, the model coupling allowing the determination of the iron losses has been applied to a 18-slots/16-poles radial-flux interior PMSM having a double-layer concentrated winding for automobile application (i.e., 16 kW @ 1,000 rpm), mainly for EVs/HEVs/FCVs.The comparisons between the model coupling results and the FE simulations have confirmed the validity of this original approach.The influence of the MEC discretization on the iron loss calculation and the electromagnetic performances has been analyzed.The gain in computation time obtained with this original approach has been given.

General Assumptions
The generalized nonlinear adaptive MEC is based on the following simplifying assumptions: • The semi-analytical model is supposed 2-D/quasi 3-D; • The end-effects are neglected; • The electrical conductivities of the materials (e.g., the PMs, the copper, and the iron) are assumed to be null; • The magnetic materials are considered as isotropic; • The mechanical stress on the nonlinear B(H) curve is ignored; • The hysteresis effects are ignored.
However, in the special case of electrical machines, the generalized nonlinear adaptive MEC accounts for: • The rotating/linear machines; • The radial-/axial-/transverse-flux machines; • The integral-/fractional-slot machines; • The internal/external rotor topologies; • The saturation, slotting and curvature effects; • The (non-)overlapping winding distribution with one or more layers; • The leakage flux; • The cross-saturation effect between d-q axes.

Automatic Mesh
In [45], Dubas et al. present a 2-D/quasi 3-D generalized nonlinear adaptive MEC using a mesh-based formulation.This semi-analytical model allows computating the global/integral quantities in the electromagnetic devices.For the electrical machines, it includes the automatic mesh of static/moving zones, the saturation effect, and the connection of the zones for the rotor motion (which is ensured by a new approach called "Air-gap sliding-line technic").This generalized semi-analytical model has been applied to an axial-/radial-flux interior PMSM [46][47][48] and a coaxial magnetic gear equipped with surface-mounted PMs [49].
In (δ, ζ, τ) coordinate system, Figure 1a represents the generalized discretization of electrical machines.We can distinguish two zones, i.e., a static zone and a moving zone, having a number of slices N slice in the τ-axis.Each zone can be independently modeled and the interconnecting between the two zones is performed by an "Air-gap sliding-line technic" [45].Consequently, the stator and/or rotor topology can be replaced by another type of construction.
Given the periodicity/symmetry of the electrical machine, the stator (rotor) can be divided into m (m ) zones according to the δ-xis and k (k ) zones according to the ζ-axis.It should be noted that the k (k ) zones can be constituted by a duplication of n (n ) zones according to the stator (rotor) structure.The air-gap is divided into a single zone according to δ-axis.The mesh elements having the same ).The single air-gap zone can be discretized into Ndag δ , which must be an odd discretization number to display the sliding-line in the air-gap.The discretization of the various zones is the same for each slice in the τ-axis.The mesh elements are so composed of BD blocks depending on the discretization chosen by the designer.Figure 1b shows an example of discretization for the mesh element {1, k − 1} with Nds δ 1 = 3 and where the number of BD blocks is equal to 6.The BD blocks, connected between them by the loop fluxes ψ l and giving to magnetic flux the possibility to flow in both directions, are described by a middle-point related to (except in the stator and rotor edges): • Four magnetomotive forces (MMFs) of branches (i.e., two δ-MMFs and two ζ-MMFs); • Four reluctances (i.e., two δ-reluctances and two ζ-reluctances) crossed by branch fluxes (or reluctance fluxes) φ l ; where l is the index of slice (i.e., l = 1, 2, . . ., N slice with N slice the number of slices in the τ-axis).

Matrix Formulation and Solving
By using the Maxwell's equations as well as the magnetic material equations, the generalized nonlinear adaptive MEC (where the loop fluxes ψ l are the unknowns) can be expressed, for each slice and ∀ζ rs (where ζ rs represents the mechanical position between the rotor and the stator), by: where l is the diagonal matrix of reluctances (of dimension N l × N l with N l the total number of branch fluxes, branch MMFs or reluctances), ψ l the loop fluxes vector (of dimension N l ψ × 1 with N ψ the total number of loop fluxes), χ l the topological (or incidence) matrix (of dimension N l ψ × N l ), and F l the loop MMFs vector (of dimension N l ψ × 1) defined by: in which MMF l represents the branch MMFs vector (of dimension N l × 1) that depend on the electromagnetic sources (i.e., PMs having any magnetization direction, winding fed with any waveform currents [50]).These matrixes and vectors are dependent upon the discretization of the various zones.Their dimensions are invariable with ζ rs [45].Knowing ψ l , the magnetic flux densities vector B l (of dimension N l × 1), for each slice and ∀ζ rs , can be determined by: where S l is the vector of the corresponding reluctance surface (of dimension N l × 1) in the BD blocks, and φ l the branch fluxes vector (of dimension N l × 1).Analytically, (1) can be solved iteratively with a constant relative magnetic permeability µ r according to the nonlinear B(H) curve at each iteration (e.g., by using the fixed point interation method).The mathematical formulation of the nonlinear B(H) curve is described by the Marrocco's function [51].The flowchart of the nonlinear system solving is illustrated in Figure 2 (for each slice and ∀ζ rs ).In this figure, µ ri represents the initial relative permeability in the magnetic circuit, and N sat is the total iterations number for taking into account the saturation effect.
Flowchart of the nonlinear system solving (for each slice and ∀ζ rs ) [45].

Description
The LS model [17][18][19] is a scalar dynamic hysteresis model.It is an H(B) model relating the applied excitation field H to the average magnetic flux density B in the cross section of the material.All the dynamic effects are considered in the model.It is based on the general knowledge of B and dB/dt terms and its identification is made with the aid of a characteristic surface H meas 0 (B, dB/dt) which is built from specific experimental measurements (under a controlled triangular waveform of B and variable frequency).In order to describe the macroscopic dynamic behavior of the material, a static hysteretic model is introduced in the main formulation.Indeed, the total field H(B, dB/dt) is expressed as: where : • The static term H stat (B, history) is dependent upon the magnetic flux density level and its history.
It is described by a simple static hysteresis model.In the latter, the static field is interpreted as a sum of an anhysteretic field and an additional one attributable to the friction of the Bloch wall displacements.• The dynamic term H dyn (B, dB/dt) is calculated from the measured dynamic surface H meas dyn (B, dB/dt).Indeed, after subtracting the static contribution from the characteristic surface H meas 0 (B, dB/dt), H dyn (B, dB/dt) is formulated into polynomial functions thanks to analytical interpolations [19].Fig. 3 shows the dynamic interpolated surface H dyn (B, dB/dt) for the investigated non-oriented SiFe steel.

Modeling and Implementation in Flux TM Software
Starting from quasi-static and dynamic hysteresis measurements made under variable frequency and waveform excitations, the LS model of the SiFe grade constituting the investigated machine has been built in Matlab ® and checked against experimental data.In almost cases, the model reproduces the magnetic behavior of the material within 10 %. Figure 4 displays some examples of measured and simulated hysteresis loops by the LS model for different magnetic flux density waveforms, different induction levels and different excitation frequencies.One observes the good agreement between the simulations and the measurements.
Once the model of the studied non-oriented SiFe steels is established and validated, it was implemented in 2-D transient FE Flux TM software as a post-processing module.Indeed, the FE simulation makes possible the evaluation of the instantaneous evolution of the magnetic flux density in each element of the mesh of the structure under test, the LS model is after that applied to post process the local and total iron losses.The LS model can thus easily be applied for the estimation of iron losses in the electromagnetic devices.It has to be mentioned that in FE simulations, only the normal magnetization curve of the soft magnetic material constituting the the magnetic structure is considered (not the entire hysteresis loop).The computation of iron losses is then performed as a post processing task [see Figure 5].

Principle
In order to exploit the outputs of the established nonlinear adaptive MEC, viz., the magnetic flux density waveform in each magnetic reluctance, a new approach allowing the calculation of the iron losses by the LS model in a BD block has been developed.This approach is then generalized on all BD blocks of the semi-analytical model.
The main idea of the new appraoch is to subdivide the BD block into four regions in which we systematically find two magnetic flux densities, viz., B l ζ and B l δ .Figure 6 shows the decomposition of one BD block for the iron loss estimation.Thus, for each region (i.e., 1/4 of the BD block) and for one temporal period of simulation, the procedure of the iron loss prediction may be summarized as below: • The assessment of the direction for which B l = B l max .This amounts to evaluate the magnetic angle ϕ l for each slice [see Figure 7], viz., • The association of a new local coordinate system to the B l max direction.It should be noted that the //-axis is taken to be horizontal to B l max and the ⊥-axis is taken to be vertical to the //-axis.• For a given instantaneous magnetic flux density B l , are associated two new components of B l , i.e., B l // and B l ⊥ [see Figure 7] which can be expressed as: • The last step consists of estimating the iron losses in both //-and ⊥-axis at every time step (or every ζ rs ) and for each slice.The new components B l // and B l ⊥ are then used as inputs of the LS model to calculate the two components of magnetic field, i.e., H l // and H l ⊥ , and deduce the corresponding hysteresis loops.The local iron losses (i.e., for each region) are determined by the enclosed area of the modeled dynamic hysteresis loop and added to give the instantaneous total iron losses in the BD block and then in the nonlinear adaptive MEC (i.e., all the BD blocks) [see Figure 7].It should be noted that the choice of the (//, ⊥, τ) coordinate system is not arbitrary.Indeed, it allows a direct application of the scalar dynamic LS model in the case of a uniaxial magnetic flux density.In addition, it is also interesting to note that the knowledge of the volume of each region is necessary for the passage from the iron losses per cycle and unit mass to the total iron losses in [W].

Description of the Machine
The coupling between the generalized nonlinear adaptive MEC and the LS model has been applied to a radial-flux interior PMSM with 18-slots/16-poles (i.e., fractional-slot number) [see Figure 8].This electrical machine is used for automotive application (i.e., 16 kW @ 1,000 rpm), mainly for EVs/HEVs/FCVs.
The topology of the machine is represented in the Figure 8a.The stator has a double-layer concentrated winding distribution (viz., the non-overlapping winding with all teeth wound), supplied by sinusoidal current waveform with a maximum amplitude of I smax .The 3-phase windings are star-connected.The winding schematic, having N stp serie turns per phase, is depicted in Figure 8a.Due to the periodicity of the winding, the electrical machine can be reduced to its half (i.e., 9-slots/8-poles).The PMs, inserted in the rotor core and radially magnetized, are isotropic and have a linear demagnetization characteristic with a maximum operating temperature T m f equal to 100 °C.The rotor bridge must be as thin as possible to avoid flux leakage.Its thickness should be obtained by optimal design to ensure best electromagnetic performances and mechanical robustness [48].The geometrical and physical parameters of the studied radial-flux interior PMSM are represented in Figure 8b and given in Table 1.

Mesh Elements and Discretization
Based on the notions of the generalized nonlinear adaptive MEC [see Section 2] and using (r, Θ, z) as a replacement to (δ, ζ, τ), a 2-D MEC of the machine (with N slice = 1) has been developed in polar coordinate system (r, Θ).Table 1.Geometrical and Physical Parameters of the Radial-Flux Interior PMSM.

Discretization Minimal Average
Nds r [1 3.9 6.7 The discretization of the static and moving zones, for the coupling between the 2-D nonlinear adaptive MEC and the LS model for the iron loss determination, is given in Table 2.The air-gap zone is discretized into Ndag r = 3.Consequently, the nonlinear system (1) is composed of: • N l ψ = 380 loop fluxes and N l = 1, 606 branch fluxes for the minimal discretization; • N l ψ = 956 loop fluxes and N l = 4, 036 branch fluxes for the average discretization.It can be solved for each rotation step and each iteration at saturation.The discretization N pts corresponding to ζ rs → Θ rs for the rotor motion is equal to 100.The total number of iterations for magnetic convergence N sat was imposed on 20 with an initial relative permeability in the magnetic circuit µ ri = 6,000.

Introduction
In order to conclude on the adequacy of the established model coupling, the iron losses of the studied electrical machine have been evaluated for the no-load (i.e., 0 A @ 1,000 rpm) and load (i.e., I smax @ 1,000 rpm) operation.The iron losses generated in the rotor zone being low, we present here only those of the stator.The influence of the discretization (i.e., minimal and average) on the iron losses and electromagnetic performances has been analyzed.The obtained results have been compared with those made retrospectively in 2-D transient FE Flux TM for the same number of rotation step (viz., N pts = 100).
The personal computer used for this comparison has the following characteristics: Intel(R) Core(TM) i7-2600 CPU@3.4GHz (with 2 processors) RAM 4 Go 64 bits.

Results Discussion
Table 3 summarizes the results of the total iron losses in the stator computed at no-load/load for the average discretization.The results reveal that the proposed approach is able to achieve a fine prediction of the iron losses from the magnetic flux density waveforms of the 2-D nonlinear adaptive MEC. Figure 10  In order to analyze the influence of the discretization on the iron losses and the electromagnetic performances, a stator tooth pitch was defined as demonstrator [see Figure 9a].The latter, subdivided into nine mesh elements, is representative of the different evolutions of the magnetic flux density in the various regions of the machine (i.e., tooth, yoke, tooth-tips...).The Appendix illustrates some  examples of mesh elements with a minimal discretization and the iron loss computation from their magnetic flux densities using the LS hysteresis model.Tables 4 and 5 summarize the comparisons of the iron loss computations in the stator under a tooth-pitch for both discretizations.It can be noted that the minimal discretization implies significant average discrepancies.These average discrepancies up to ≈ 15% are observed at load and up to ≈ 28% at no-load [see Table 4].The examination of the magnetic flux density waveforms versus Θ rs in Figure 11 allows the explanation of these discrepancies due to this coarse mesh size.Indeed, focusing, for example, on the case of the tooth (the mesh element no.6 of Figure 9a ), at load, the rand Θ-components of the magnetic flux density at load provided by the 2-D nonlinear adaptive MEC (used to perform the iron loss calculation in the stator) are different from those resulting from the FE simulation of the electrical machine.It should be noted the data acquired in the Flux TM environment are associated to the point sensors (spatial magnetic flux density on a given point) whose the geometrical coordinates correspond to the middles of the magnetic reluctances.In terms of electromagnetic performances, the average electromagnetic torque Γ em for the minimal discretization is equal to 50 Nm, much different in comparison with that of the FE simulations (viz., 150 Nm).
For the average discretization, the results show the good evaluation of the local iron losses in the nine mesh elements from the magnetic flux density waveforms of their associated BD blocks.Notice that the average total iron losses in the stator under a tooth-pitch are computed within ≈ 5 % in comparison with the calculations obtained by the FE simulations [see Tables 5].In Figure 12, always for the same mesh element no.6 (i.e., the stator tooth), it can be observed that the r-component of the magnetic flux density at load provided by the 2-D nonlinear adaptive MEC (used to perform the iron loss calculation in the stator) is quite similar to that of the FE simulation.However, it is interesting to note a light difference between the two Θ-components of the magnetic flux density.The iron loss calculation was not impacted by the Θ-component due to the lower values of magnetic flux density in comparison to the r-component and to the change of coordinate system (i.e., // et ⊥) adopted during the calculation.The electromagnetic torque with the average discretization is much the same as the reference one given by the FE simulations.

Conclusion
In this work, a coupling between a 2-D/quasi 3-D generalized nonlinear adaptive MEC (using a mesh-based formulation) and the LS hysteresis model has been established.This original approach allows the computation of the magnetic losses in the electromagnetic devices from the magnetic flux density waveforms of the MEC.
In order to evaluate the reliability of this approach, the developped model coupling has been applied to a 18-slots/16-poles radial-flux interior PMSM having a double-layer concentrated winding, dedicated to automotive applications (i.e., 16 kW @ 1,000 rpm), mainly for EVs/HEVs/FCVs.Two discretizations (viz., minimal and average) of the 2-D nonlinear adaptive MEC of the electrical machine have been adopted.The influence of the discretization on the iron loss evaluations and the electromagnetic performances has beed studied.The results of the model coupling have been compared with those made retrospectively in the 2-D transient FE Flux TM .In conclusion, this approach allows reducing computation time in comparison with the EF simulations, while providing quite similar results with a sufficient discretization of the mesh elements of the 2-D nonlinear adaptive MEC.

Figure 1 .
Figure 1.Generalized discretization of electrical machines [45]: (a) Automatic mesh of electrical machines for each slice in the τ-axis, and (b) Discretization of a mesh element (e.g., for {1, k − 1} in the stator).

Figure 3 .
Figure 3.The dynamic interpolated surface H dyn (B, dB/dt) of the investigated non-oriented SiFe steel.

Figure 7 .
Figure 7. Change of the coordinate system for the magnetic flux densities & principle of the iron loss calculation.

Figure 9 .
Figure 9. Mesh elements of the radial-flux interior PMSM in (r, Θ, z) coordinate system: (a) Static zone under a tooth-pitch, and (b) Moving zone under a pole-pitch.
displays the computation time for the 2-D nonlinear adaptive MEC which is divided by 3/2 in comparison with the 2-D transient FE Flux TM .

Figure 11 .
Figure 11.Magnetic flux density waveforms at load for the minimal discretization.

Figure 12 .
Figure 12.Magnetic flux density waveforms at load for the average discretization.

Table 2 .
Discretization Studied for the Model Coupling.

Table 3 .
Total Stator Iron Losses for an Average Discretization.

Table 4 .
Iron Losses in the Stator under a Tooth-Pitch for the Minimal Discretization.
* see Figure9afor the numerotation of the mesh elements.

Table 5 .
Iron Losses in the Stator under a Tooth-Pitch for the Average Discretization.
* see Figure9afor the numerotation of the mesh elements.