Fast and accurate methodologies for the eddy current losses analysis in axially laminated rotor of synchronous reluctance machines under different loads

The authors present two fast and accurate methodologies for the computation of eddy current losses in the axially laminated rotor of a synchronous reluctance machine. The methodologies are based on different combinations of the finite element method in time and frequency domains with 2D and 3D formulations. First, a comparative study of the 2D and 3D formulations for loss calculation is presented, considering various load angles of the machine to illustrate the problem of eddy ‐ current losses in this type of machines and its dependence on the load angle. The influence of the iron saturation on the loss calculation is also evaluated in these computations. A novel correction factor based on the computations at two load angles is proposed to convert the losses computed from a 2D model to match those computed from a 3D model. For the sake of generality, investigations are also conducted for various thicknesses of the lamination layers and different machine lengths, and an analytical method to describe the dependency of eddy ‐ current losses on the load angle of the machine is introduced. Moreover, a simplified method is proposed for modelling eddy currents in the frequency domain and calculating losses in an axially laminated structure based solely on the results of a magnetostatic solution. The results obtained by the simplified model demonstrate excellent agreement with the full 3D magneto ‐ dynamic simulation. Overall, the findings contribute to understanding and accurately characterising the eddy current losses in axially laminated rotors, offering potential insights for designing and optimising axially laminated synchronous reluctance electric machines.


| INTRODUCTION
High-speed synchronous reluctance machines with axially laminated rotors have received limited attention in the literature, and a well-established methodology for their design is yet to be developed [1].During the design stage of the machine, it is crucial to accurately and efficiently calculate the losses to ensure the correct sizing and effectiveness of the cooling system.While several widely used models for iron loss calculation, such as the Steinmetz or Bertotti model [2], exist, they are primarily suitable for machines with transversely laminated rotors and relatively thin lamina at low frequencies.In the case of synchronous reluctance machines with axially laminated anisotropic (ALA) rotors, applying these empirical equations to predict losses leads to inaccurate calculations, especially at high frequencies, as will be demonstrated with the numerical simulations presented in this paper.
Previous studies that explored the issue of eddy currents in ALA rotors [3][4][5] primarily focused on the 2D formulation and provided analytical models applicable solely to the studied This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.cases.In ref. [4], the authors emphasised the necessity of conducting a comprehensive 3D analysis to accurately calculate eddy current losses in rotors with this particular configuration.However, due to the substantial computational resources' requirements associated with such an analysis, it has not been undertaken thus far.The authors in ref. [7] proposed introducing a correction factor between the 2D and 3D calculations of eddy current losses, thereby avoiding the computationally expensive 3D analysis at every iteration step of the design procedure.However, the technique described in ref. [7] is not applicable to axially laminated rotors as it assumes uniform saturation levels and skin depths for permanent magnets.In ALA machines, the position of the rotor (mechanical load angle) directly influences the saturation conditions and the depth of the field penetration into the ferromagnetic regions of the rotor, which needs to be accounted for when seeking accurate computations.To our knowledge, no method in the literature combines different computation formulations to alleviate the computation time while keeping the same level of accuracy as the full 3D magneto-dynamic results.This paper presents two different methodologies for accurately calculating the eddy current losses in ALA synchronous machine-type, thus obviating the need for a complete 3D magneto-dynamic analysis (time-stepping).The accuracy and computational efficiency of the proposed methods are demonstrated by comparing their results to theses obtained from complete 3D magneto-dynamic computations as well as by comparing the proposed methods to each other.
The paper is organised as follows: in section 2, the machine's specification is given, and different features and concepts related to the eddy current losses are presented and illustrated with results from 2D and 3D computations.In section 3, the proposed methodology is introduced; it consists of two alternative workflows combined with an approximation function for eddy-current losses calculation for the entire loading range of the machine.In section 4, the results from both workflow alternatives are presented and compared to each other.Additional 3D computations are also made to validate the results at several machine operation points.The validation is performed with computation results that have not been used in the fitting process of the approximation function.Section 5 discusses the results from the accuracy and computational efficiency points of view.

| Machine specification
The investigations in this work are focused, without losing the generality, on a 290-kW high-speed synchronous reluctance machine with an ALA rotor, the configuration of which is shown in Figure 1.It should be noted that in all the computations, the machine was supplied with a sinusoidal current source.
The main parameters of the electric machine are given in refs.[8,9] and are reported in Table 1.

| Machine description
As shown in Figure 1, this type of rotor presents a path for the eddy current to flow in the axial direction, making the eddy current losses problematic, particularly for high-speed applications.It is worth highlighting that the slot harmonic is the primary harmonic component responsible for inducing eddy current losses in the rotor as depicted in Figure 2.These results are obtained from a complete 3D magneto-dynamic computation, the details of it are described in the methodology section.
Understanding the dominance of this harmonic component is crucial for accurately predicting and mitigating eddycurrent losses in the rotor.As it can be seen from Figure 2, the highest amplitude of the current density is observed at 18,900 Hz, which exactly corresponds to the frequency of the slot harmonic expressed as the following equation: where Q s is the number of stator slots, f 1 is the supply voltage frequency, and p is the number of pole pairs.

| METHODOLOGY FOR THE EDDY CURRENT ANALYSIS
Figure 3 shows the calculated values of rotor losses obtained using complete magneto-dynamic 3D and 2D modelling.They were obtained for the steady state operation conditions where the maximum losses occur in the no-load condition and at 2 points in the operating range (20°and 30°).It should be noted that the load angle in this work refers to the mechanical load angle.
As can be seen, the 2D approach gives a significant discrepancy compared to the 3D approach.However, the calculation in 3D takes about 6 h of computing time for a single loading operation, whereas the 2D computation took no more than 22 min, including the magnetostatic computation to initialise the time-stepping and avoiding unnecessary computational transient.Calculating losses over the entire loading range will be unnecessarily long and computationally intensive.The main objective of this study is to propose two alternatives for fast and least expensive calculation of losses over the entire load range.
We propose two alternatives to calculate the eddy current losses in the ALA rotor of synchronous reluctance machines for the entire load angle range, which are illustrated in Figure 4.The first alternative consists of carrying out a complete magneto-dynamic 3D computation at two loading points, from which a curve fitting is extracted for the rest of the loading range.The second alternative is similar to the first one, but instead of complete 3D magneto-dynamic computations we propose to make a magnetostatic 2D computation at two loading points and, from each of them, extract the slot harmonic in the airgap through Fast Fourier Transofrmation (FFT), which is then used in the harmonic 3D analysis of the rotor only.The results of these two loading points are then used to fit the losses curve and extract them for all the other loading conditions.In both alternatives, knowledge of only two loading points is required to approximate the entire range of loads, corresponding to the minimum and maximum losses.These loading points can be obtained through 2D computations as they follow the same trend as in the 3D.
The following sections will explain the details of the computations and the approximation and show their efficiency and accuracy.

| FEM formulations for complete 2D and 3D modeling
For 2D analysis, the formulation of magnetic vector potential (A-formulation) was used.In this case, fine mesh full-timestepping analysis with the initial magnetostatic step takes about 22 min.However, a complete 3D analysis is associated with mapping problems between a stationary and moving mesh.To avoid convergence problems, it is thus necessary to use a mixed vector-scalar formulation (A-φ formulation).Complete time-stepping analysis in 3D with initialisation takes about 6 h for one period.Figure 5 shows the results of the magnetic flux density distribution for 2D and 3D magnetostatic solutions at the no-load point.

| Formulation of the eddy current
According to ref. [6], using Gauge fixing, the magnetic diffusion equation in terms of magnetic vector potential is given by the following equation:

F I G U R E 4
The two proposed alternatives to compute losses in the axially laminated anisotropic (ALA) rotor.
SITNIKOV ET AL. -3 where σ is the material's conductivity, Α the magnetic vector potential, J s the source current density, and μ the magnetic permeability of the material.Once ( 2) is numerically solved, the eddy current losses can be calculated as follows [4]: where the integration is made over the volume of the conducting region of the rotor (V).

| Approximation of losses
Based on the simulation data, the following approximation equation for the eddy current losses at any loading angle is proposed: where P max and P min are maximum and minimum eddy-current losses, respectively, sech is the hyperbolic secant function, p the number of pole pairs, α the mechanical load angle in radians and α m−g the transition angle from motor to generator mode (positive torque changes to negative) in radians.

| 3D harmonic rotor model (3D HRM)
A novel 3D harmonic rotor model (3D HRM) is introduced to streamline and expedite the modelling and calculation of eddy current losses, relying solely on the outcomes of the magnetostatic analysis.
The proposed 3D HRM operates on the following assumptions: 1.It considers that the slot harmonic exclusively causes the eddy current losses.2. The magnetic permeability is assumed to be constant, the value of which is locally set according to the results of the 2D magnetostatic analysis (frozen permeability method).
3. The influence of the end winding is disregarded in the calculations.
By employing these simplifications, the proposed 3D HRM achieves efficiency and computational expedience while still providing valuable insights into eddy current losses in the rotor.Furthermore, excluding the end winding effect contributes to the overall simplification of the model, ensuring a more straightforward and efficient analysis process.
The concept of a simplified model is presented in Figure 6.
The proposed simplified model offers several key advantages.Firstly, it enables the simulation of only the rotor and a portion of the air gap, eliminating the need to model the stator and winding.This approach improves the mesh quality and facilitates more accurate analyses explicitly focused on the rotor region.Another significant advantage of the HRM is its utilisation of a frequency domain solver without time-stepping analysis, significantly reducing computing resources.This streamlined approach simplifies the modelling process, reducing complexity and computational overhead.The main physics components employed in the simplified model are as follows: 1. Airgap: The airgap is modelled, and the source is imposed as a magnetic flux density boundary condition, representing the excitation applied to the system.The source itself is extracted from a 2D magnetostatic computation of the whole machine, followed by an FFT of the magnetic flux density in the airgap of the machine.

| FFT analysis of the airgap harmonics
3D HRM includes several essential but, at the same time, understandable steps, which will be described in this paragraph.This model consists of the 2D magnetostatic analysis of two load points from which slot harmonic components of magnetic flux density are extracted and used as a source in the airgap.Permeability is also calculated in 2D magnetostatic cases, extracted as a function of coordinates in XY-planes and used for steel rotor domains in HRM.The general workflow with the 3D HRM is presented in Figure 7 and main equations are as follows.
Stator slot harmonics from FFT analysis can be expressed as the following equation: Stator slot harmonics in rotor coordinate frame: Stator slot harmonics in the frequency domain (exponential form):

| Mesh configuration
It is essential to acknowledge that the unique characteristics of a high-speed electric machine without traditional lamination of layers pose specific challenges when dealing with field penetration depth.To ensure accurate analysis, a dedicated fine enough mesh is necessary.However, balancing computational resources and result accuracy becomes crucial while refining the mesh.Figure 8 illustrates the default mesh generated by the software, utilising the machine's geometry and selected physics.
One of the criteria for assessing the quality of the mesh is the skewness, which is expressed as follows [10]: where θ is the angle over a vertex (2D) or edge (3D) in the element, θ e is the edge or vertex angle in an ideal element and the maximum is taken over all vertices (2D) or edges (3D) of the element.In the 2D case, the average mesh quality is 0.76 (with a maximum value of 1) for a mesh of 8500 elements.
Increasing the number of elements to refine the mesh can improve the accuracy in computing the eddy currents; however, this will lead to a significant increase in the computing resources.The mesh configuration we propose can strike a balance between accuracy in computing eddy currents and computational resources, such a mesh is shown in Figure 9.
This mesh consists of mapped elements inside the layers and triangular elements on all other parts.In turn, the mapped elements have an exponential symmetric distribution over the thickness and length of the layer.The average quality of the mesh was 0.82 using 19,275 elements.Figure 9 shows the computed eddy current density distribution using the proposed mesh.
For the 3D analysis, we maintain the mesh quality by employing a swept mesh technique, which extrudes the original mesh along the length of the machine and generates a 3D mesh that accurately represents the electromagnetic phenomena while preserving the computational efficiency.

| Mesh
The first results are related to the mesh quality used in 2D simulations.Despite the high-quality mesh and low computational time (70 s), it has been observed that the distribution of eddy current density is not correctly computed (Figure 10).This can be seen as a uniform distribution of eddy-current density between the rotor layers corresponds to the stator slot opening.However, we should expect that the distribution of eddy current density exhibits low values in the ferromagnetic regions, where the field penetration depth is small, and more noticeable values in the non-magnetic regions of the rotor.
The inaccuracies in the eddy current density distribution indicate the limitations of the default mesh in capturing the complex electromagnetic phenomena at play.It emphasises the necessity for a more refined and dedicated meshing approach to achieve accurate and meaningful results.
Indeed, the current density distribution obtained with the proposed mesh configuration presents a physically substantiated picture (Figure 11).This accurate representation of the current density supports capturing the actual distribution of eddy current behaviour within the axially laminated rotor.To assess the need for a refined mesh, we conducted a sensitivity analysis, as reported in Tables 2 and 3, for rotor layers of 1 mm thickness.
F I G U R E 9 Proposed mesh with a mix of mapped and triangular elements.
F I G U R E 1 0 Eddy-currents density using a default mesh (greyelectrical steel, blue-non-magnetic steel Inconel).

F I G U R E 1 1
Eddy-currents density using the proposed mesh (greyelectrical steel, blue-non-magnetic steel Inconel).

Elements per layer
Eddy-losses, W Time, s Table 2 compares the eddy current losses for different numbers of elements across the thickness of the layer with 40 elements per arc length, as well as the time spent on modelling the whole rotation period.
Table 3 compares the eddy current losses for different numbers of elements across the arc length of the layer with 4 elements in the thickness and the time spent modelling the whole rotation period.
From Tables 2 and 3, it can be seen that with three elements in thickness and 40 elements in length, the calculation process begins to stabilise, and the results differ little with further refinement; however, the calculation time increases significantly.Thus, 3 elements along the layer thickness and 40 elements along the length are necessary and sufficient for correctly calculating the rotor eddy current losses.
It is important to highlight that COMSOL swiftly generates the default mesh within 1-2 s, whereas the proposed mesh requires manual configuration.The setup time for the mesh involves manually identifying and marking all necessary rotor layers, as automation possibilities are limited.However, this manual process is accomplished within a concise timeframe of no more than 5 min.Once the rotor mesh is established, COMSOL can automatically generate the remaining components.Consequently, the time invested in mesh setup can be integrated into the overall model computation time, with the advantage that the mesh configuration only needs to be performed once, without necessitating further adjustments for subsequent calculations.

| 2D losses calculation
A complete magnetostatic and magneto-dynamic (time-stepping) study was done in 2D. Figure 12 shows, for a layer thickness of 1 mm, the eddy current losses as a function of the load angle; the per unit torque is also shown in the figure.

| Losses dependency on the lamination thickness
Several rotor lamination thicknesses have also been tested at the same saturation level, and they all have the same trend, which can be described in Figure 12 and exact values could be obtained by multiplying the pu values from the figure with values from Table 4.

| Results of approximation
Approximation curves were obtained based on the method described above using equations (6). Figure 13 shows the dependencies of approximation and losses calculated by 3D magnetodynamic FEM.

| Simulation results of the 3D harmonic rotor model
Figure 14 shows the distribution of eddy current losses at the load angle of 45°(where the maximum losses are observed).
Knowing the distribution function of losses depending on the load angle, we need to conduct studies only for two rotor positions for a simplified 3D model-no-load and transition from motor to generator mode-and then we can get all the missing losses using the approximation function.
Figure 15 compares the results of a complete 3D calculation using an approximation and a calculation using a simplified 3D formulation also using an approximation.SITNIKOV ET AL.
- A 3D HRM makes obtaining results quite close to 3D magneto-dynamic modelling possible.At the same time, for a complete 3D calculation with a good mesh, it takes 6 h to calculate one period, while for a simplified 3D model, the calculation takes 1 min.
Table 5 compares the maximum rotor eddy current losses per pole of different machine lengths when calculated using 3D HRM and full 3D time-stepping models.
Table 6 compares the maximum and minimum losses per pole with fixed lengths for different rotor layer thicknesses.
As can be seen from Table 6, there is a slight discrepancy at 2 mm thickness.This can be explained by the fact that, in this case, due to the alternation of layers, the top layer turns out to be ferromagnetic and creates an additional reaction field in the HRM, which is not compensated (Figure 16).
Table 7 compares the maximum rotor eddy current losses per pole of different machine lengths with 5 mm lamina thickness when calculated using 3D HRM and full 3D timestepping models.
F I G U R E 1 5 Normalised losses comparison for complete 3D timestepping analysis and using simplified 3D method (2464 W corresponds to 1 for eddy current losses).As can be seen, although the maximum computational error is slightly higher for Alternative 2, the gain in computational time is the main advantage of this alternative.Also, for this case, the iterative process HRM can be implemented to take into account the reaction field and reduce the computational error, which will be done in future works.
Presently, a prototype for validating the losses computed using 3D FEM is not available.However, the mere existence of a prototype does not ensure accurate losses calculation, as they cannot be directly obtained but rely on indirect methodologies, such as temperature measurement.In order to cross-validate the results obtained using 3D FEM, a simulation of a comparable 2-pole machine with an ALA rotor, as described in ref. [3], was conducted.According to the article, with an output power of 60 kW at 48 krpm, the losses would be in the range of 228-464 W (depending on the winding configuration).Our FEM model yields a value of 283 W for an output of 71 kW, aligning with the specified range.Consequently, it can be asserted that the model functions accurately.

| Approximation
As shown in graphs 13 and 15, the algorithm gives a good approximation for the eddy current losses approximation.For a sufficiently accurate approximation of eddy current losses, it is sufficient to perform a 2D magnetostatic study for two load angles, the no-load (0 deg.) and the transition between motor and generator mode (45 deg.).
It is important to note that based on the model data, it can be seen that in 3D, the losses have the same distribution function depending on the load angle as in the 2D case; however, there is a difference between the calculations.In the no-load case, the difference between the 2D and 3D is 27%, and at the maximum load the difference is virtually zero.Knowing this difference allows us to introduce a dynamic correction factor to convert losses from 2D to 3D.

| 3D harmonic rotor model
The simplified model achieves a comprehensive system representation by incorporating these physics components while significantly reducing the computational effort required.This approach enables efficient and accurate analyses, focussing specifically on the relevant aspects of the rotor and air gap without the need to simulate the entire system.
As can be seen from the above comparisons, the simplified model makes it possible to obtain close results while significantly reducing computational costs.A slight discrepancy can be observed at a thickness of 2 mm, which can be explained by the fact that due to the symmetrical laying of the layers, the electrical steel layer became the last on the rotor, that is, it plays the role of a pole shoe.Table 2 compares the maximum and minimum loss values since it was shown earlier that knowing these two values makes it possible to describe all other values accurately.Thus, if we guarantee the accuracy for extreme values, we can guarantee the accuracy for all other cases.

| CONCLUSION
Our study has revealed a significant correlation between the eddy current losses of an axially laminated rotor and the load angle.We have proposed an analytical expression that provides a reliable approximation of these losses, requiring only knowledge of the maximum and minimum loss values.Our approximation has been validated for various rotor lamination axial thicknesses, demonstrating promising results with low error.An essential aspect of our investigation was the emphasis on 3D analysis of the electric machine.The comparison between 2D and 3D models has highlighted the substantial differences in loss calculations, with variations of up to 30% observed in no-load conditions.Furthermore, we have demonstrated that our proposed formula accurately describes the losses obtained from 3D analysis.To address the computational costs of 3D modelling, we have developed a 3D harmonic rotor model utilising a magneto-harmonic solver in the frequency domain.This streamlined model allows for efficient and accurate calculation of the two necessary points for the approximation, enabling the estimation of losses throughout the entire load range.The harmonic model relies on data from the 2D magnetostatic solver and requires only the spatial FFT analysis of the flux density in the airgap.These findings provide valuable insights for designing and optimising synchronous reluctance electric machines, particularly in highperformance applications that demand prompt and accurate analysis.
Moving forward, our future work will explore the challenges of axially laminated rotor insulation and the cosimulation modelling of the machine with a power converter.These endeavours aim to enhance further the understanding and performance of such machines in practical applications.

©
2024 The Author(s).IET Electric Power Applications published by John Wiley & Sons Ltd on behalf of The Institution of Engineering and Technology.IET Electr.Power Appl.2024;1-10.wileyonlinelibrary.com/journal/elp2 -

F I G U R E 1 F I G U R E 2 2 -
Axially laminated rotor construction (grey-electrical steel, blue-non-magnetic steel Inconel).The shaft of the machine is not shown.T A B L E 1 Main machine specification.Eddy current density frequency spectrum at a given point in the rotor, computed from a complete 3D magneto-dynamic FE analysis.SITNIKOV ET AL.

F
I G U R E 3 2D and 3D FEM computed rotor losses at different load angle showing the maximum, minimum and intermediate values (2464 W correspond to 1 in pu).

2 .
Passive conductor condition: This condition is applied to all steel rotor layers, ensuring the closure of the eddy current lines within all rotor (I c = 0).3. Magnetic flux conservation: Magnetic scalar potential (φformulation) for airgap domain.4. Zero potential at the outer boundary point: The outer boundary point of the model is set to zero potential, establishing a reference point for the analysis.

F I G U R E 5 F I G U R E 6 4 -
Magnetic flux density distribution computed with 2D and 3D magnetostatic analysis at no-load.Concept of the 3D HRM.SITNIKOV ET AL.

F I G U R E 7 5 17518679, 0 ,
3D harmonic rotor model workflow.F I G U R E 8 Default mesh generated by the software with additional conditions.SITNIKOV ET AL. -Downloaded from https://ietresearch.onlinelibrary.wiley.com/doi/10.1049/elp2.12450 by Aalto University, Wiley Online Library on [04/06/2024].See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions)on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons LicenseAll simulation results presented in this section were obtained using COMSOL Multiphysics software.The work was carried out on a PC with the following specifications: � Processor -12 th Gen Intel(R) Core (TM) i5-12,600K 3.70 GHz � Installed RAM -32 GB � System type -64-bit operating system � Windows -Windows 10 Education

3
Mesh sensitivity on the number of elements in arc length.

F I G U R E 1 2F I G U R E 1 3
Normalised eddy current losses and torque as a function of the load angle (89 Nm corresponds to 1 pu for the torque, and 2464 W corresponds to 1 pu for the eddy current losses).T A B L E 4Maximum torque and maximum losses as a layer thickness function.Losses approximation for 3D case (2464 W corresponds to 1 for eddy current losses).

T A B L E 5 6 T A B L E 6 4 F I G U R E 1 6
Comparison of losses for the complete and simplified model as a function of the machine length at the maximum eddy current losses.Comparison of losses for the complete and simplified model as a function of the layer thickness at the maximum eddy current losses.Volumetric loss density from the 3D harmonic rotor model.Volumetric loss density for a lamina thickness of 2 mm.T A B L E 7 Comparison of losses for the complete and harmonic rotor model (HRM) as a function of the machine rotational speed at the maximum eddy current losses.

Table 8
compares computational error and time for the two alternatives used.