Node Mapping Criterion for Highly Saturated Interior PMSMs Using Magnetic Reluctance Network

Interior Permanent Magnet Synchronous Machine (IPMSM) are high torque density machines that usually work under heavy load conditions, becoming magnetically saturated. To obtain properly their performance, this paper presents a node mapping criterion that ensure accurate results when calculating the performance of a highly saturated IPMSM via a novel magnetic reluctance network approach. For this purpose, a Magnetic Circuit Model (MCM) with variable discretization levels for the different geometrical domains is developed. The proposed MCM caters to V-shaped IPMSMs with variable magnet depth and angle between magnets. Its structure allows static and dynamic time stepping simulations to be performed by taking into account complex phenomena such as magnetic saturation, cross-coupling saturation effect and stator slotting effect. The results of the proposed model are compared to those obtained by Finite Element Method (FEM) for a number of IPMSMs obtaining excellent results. Finally, its accuracy is validated comparing the calculated performance with experimental results on a real prototype.


Introduction
The demand for Permanent Magnet Synchronous Machines (PMSM) is rapidly increasing in high-performance applications, such as electric vehicles [1][2][3][4], due to their high power density.In particular, the configuration of buried magnets inside the rotor is becoming very popular, because of the additional torque made available to saliency, their wide constant power speed region, and their high demagnetization withstand capability, among others [5,6].
The design process of a PMSM frequently involves the aid of software based on Finite Elements Methods (FEM) [7][8][9][10][11][12][13].Although accuracy of the results is very high, it requires a high computational cost, together with an elevated amount of time to define the problem.This makes FEM more suitable for validation purposes rather than for preliminary machine design by iterative process.Consequently, the use of analytical design tools that rely on magnetic circuits instead of FEM is becoming widespread.
Different authors have presented simple magnetic circuits for different PMSM topologies [14][15][16].The results are acceptable as a first estimation of the machine performance, but they are not comprehensive enough if a transient analysis or a deeper study is required.For this purpose, complex magnetic circuit models based on magnetic reluctance network, known as Magnetic Circuit Model (MCM), have been proposed [17].The methodology is very similar to simple magnetic equivalent circuits; the main difference lies in the larger number of elements that the machine's geometry is discretized.
In the literature, different MCMs have been proposed [18][19][20][21][22][23].However, a clear meshing criterion that guarantees accurate results regardless of the size or geometry of the machine, i.e., a node mapping criterion for different machine regions, is not provided.In general, the reluctance element's distribution is set according to the main flux paths [20,24].However, in machines that work at heavy load conditions, the flux paths are in most cases unpredictable, especially in IPMSMs, where the complex geometry of the rotor makes difficult their analytical modelling [8,25,26].In addition, the fact that the magnetic flux can only pass through an element in a unique defined direction makes it advisable to establish a general node mapping criteria for IPMSMs.
In this paper, a general node mapping criterion for IPMSMs of any geometry and size is presented.To that end, MCM with variable discretization levels for the different geometric domains is developed.The proposed MCM models V-shaped IPMSMs, with variable magnet depth and angle between magnets.A suitable MCM structure composed of generic cells, named nodal elements, as modelling element for any part of the machine is proposed.The nodal elements contain geometric and electromagnetic information of the modelled physical domain.The model contemplates rotor motion, allowing dynamic analysis, i.e., time stepping simulations.The MCM takes into account the effects of stator slotting, the airgap magnetomotive force (MMF) waveform due to armature current, and cross-magnetizing saturation effects due to an improved magnetic flux path definition.This, in turn, leads to more accurate results than other simpler magnetic circuit models shown in the bibliography.
This paper is organized as follows.First, the MCM developed for IPMSMs is described.Next, the magnetic phase flux linkage at heavy load conditions for different V-shaped IPMSMs is calculated.Then, an analysis of the accuracy of the results from using different node mappings in the MCM is conducted, and the results are compared to those obtained by the FEM software.This analysis provides a node mapping criterion to ensure sufficient accuracy when calculating the performance of an IPMSM with a magnetic reluctance network.Finally, once the node mapping criterion is described, the presented methodology is validated by comparing the calculated performance to tests on a real traction IPMSM working under heavy load conditions.

Definition of the Magnetic Circuit Model
In this section, the MCM developed for IPMSM is described.The MCM allows the number of nodes the machine is discretized into to be selected.The proposed MCM is divided into three regions: stator, rotor and airgap.A generic nodal element, shown in Figure 1, is proposed.The node is set in the center of the nodal element.Each nodal element has its own coordinates (i,j), so it may be stored in a nodal elements matrix.Appropriate row and column are assigned to each nodal element according to its location on the associated discretized geometrical domain being modelled.Information regarding the associated geometry (identifier, row, column, region, position, etc.) is stored as well.Each nodal element has four sub-elements that capture both the radial and orthoradial magnetic fluxes and, therefore, possible cross-coupling effects [27,28].The sub-element geometrical domain can be rectangular, trapezoidal or even triangular.Sub-element information that is held is fundamental for the solving process: orientation length L, transversal section A i , permeance, permeability µ, MMF source, etc. Magnetic permeance, P i,orientation , is employed for convenience in the solving process; it is the inverse of magnetic reluctance, i,orientation , which is calculated by Equation (1) [29].
To reduce the computational effort, if periodicities exist within the machine, just part of it is modelled.The number of periodicities is obtained according the number of pole pairs and slots (Equation ( 2)).
where Q s is the number of stator slots, p the number of pole pairs and GCD stands for Greatest Common Divisor.
The geometric positioning of the nodal elements that model different regions of the machine, is established according to the existing main magnetic flux paths.In addition, it facilitates the connection between nodal elements, and therefore, the different regions of the machine can be accurately meshed.

Stator
The stator yoke, slots, slot-openings, teeth, and tooth-tip are modelled separately.The MCM parameters that define the stator node mapping are presented in Table 1.As an illustrative example, the values given in the table result in the node mapping shown in Figure 2. Presented parameters in Table 1 define the number of node columns and rows in stator regions, as can be observed in Figure 2.For open slot machines, N col,t0 is null.The armature current MMF sources are located at the north and south sub-elements of each element belonging to the teeth and are calculated by Equation (3).
where [Ws] is the winding sequence matrix that relates the teeth wound by each coil to the corresponding phase winding [30].[I UVW ] is the phase current vector.Z slot,ph,s is the number of turns per slot and layer.Finally, N w,paral,s is the number of parallel connected winding groups.

Rotor
The parameters that define the reluctance network for IPMSM rotors are presented in Table 2. Due to the existing symmetry, only the mapping for half a pole needs to be defined.To ease the comprehension, the values given in the table result in the node mapping shown in Figure 3.In Figure 4, the whole pole node mapping scheme for V-shaped IPMSMs is shown.Regarding the rotor nodal mapping, some nuances must be taken into consideration: • N row,PM also define the number of node rows in the non-magnetic material and in the q-axis magnetic bridge.

•
N col,PM : -It defines rotor yoke's node mapping: over and below PM side.These node distributions are "triangular", with N col,PM columns and as we move away from d-axis, the number of rows decreases from N col,PM to one.-It defines with other parameters as N col,bridge and N col,bridge q and one necessary extra node (for node connection purposes), the number of columns in the N row,bridge upper rotor nodal rows, belonging to an arc whose thickness is the bridge height.
It is remarkable that the same nodal mapping defined by Table 2 may be established for extreme cases of IPMSM.As an example, V-shaped in Figure 5    The magnets, the rotor yoke, the magnetic bridge and the rotor slot non-magnetic material, which is responsible for preventing a large magnetic flux leakage, are modelled by a selectable number of elements.The MMF contribution is computed by Equation ( 4) and assigned to the north and south sub-elements corresponding to permanent magnets.
h PM is the magnet height, and H c the coercive field strength.

Airgap
This region is the most significant part of the MCM and the overall precision of the model depends on its modelling.Moreover, it is the link between the stator and rotor models, thus a correct modelling of the airgap is necessary.Therefore, the airgap model is built by taking into account the relative position between the stator and the rotor.Since machine rotation is considered, airgap node mapping is varied for each time step.An airgap element is placed at the same angular position of each stator and rotor nodal element in contact with the airgap, as shown in Figure 7.It is important to note that more finely discretized stator and/or rotor models entail a more detailed airgap region node mapping.
The number of airgap row elements is determined by the parameter N row,agap .In Figure 7, the proposed airgap node mapping is visually displayed for a N row,agap equal to 2. To link the airgap with stator and rotor entities, auxiliary nodes are placed in the airgap region boundaries, as shown in Figure 7 with small yellow circles.They only have north and south sub-elements, and an infinite permeance.Therefore, it is a mathematical element that is introduced in the circuit matrix system.Its necessity is due to the inter-entities borders, where the airgap sub-elements are connected to the nearest auxiliary node (Figure 8).If they are connected simultaneously to various airgap nodes, the existence of these auxiliary nodes guarantees that each nodal element belonging to any machine entity is connected to other nodes by no more than four sub-elements.Thus, it is possible to generate the flux and permeanace matrix, according to the defined nodal element structure (Figure 1).

Solving Process
Once all the nodes are linked and the branches defined, the MCM is completely set.Owing to the non-linear magnetic behavior of the stator and rotor materials, whose B-H curve are imported from materials database and used as a lookup table, the MCM needs to be solved iteratively.Firstly, permeance and magnetic source flux matrices, [P] and [φ], are calculated for each sub-element.The matrix circuit is solved in terms of Kirchoff's Voltage Law (KVL), and the scalar magnetic potential at each node is obtained by: The magnetic flux that crosses each corresponding branch from node i to node j is calculated by: where MMF(i, j) is the addition of different MMF sources at each branch.A weighted average of the updated and previous iteration magnetic permeability is used in the following iteration step by the corresponding permeance sub-element, rewriting data in [P] and recalculating Equations ( 5) and ( 6) [20].The iteration process is finished when the following convergence criterion is met: where µ error,k is the committed mean error in sub-element permeability at iteration k.

Data Post-Processing
A key parameter for determining machine behavior is the phase magnetic flux linkage, Ψ ph [31].The magnetic flux linked by each phase winding can be computed by Equation (8).
where N sim is the number of symmetries which the problem has been divided due to the existing geometry periodicity, Z slot,ph,s is the number of turns per slot and layer, and N w,paral is the number of pole groups in parallel.The matrix [Φ t,s ] is the magnetic flux that crosses each pair of slot-teeth.This is easily obtained once the MCM is solved.Depending on whether the MCM is solved for a load operating point or at no-load condition, the phase magnetic flux linkage is denoted as Ψ load,ph and Ψ PM,ph , respectively.Using Park's transformation matrix, it is possible to work at d-q rotational reference system, whose main advantage is the fact that the different variables are time invariant.The electromagnetic torque is computed by Equation ( 9).
The phase back EMF (Equation ( 10)) and the phase voltage (Equation ( 11)) are obtained by the time derivatives of the no-load and load fluxes: where I ph is the phase current and R ph is the stator phase winding resistance.

Results Using MCM with Different Node Mapping
The described MCM is implemented in MATLAB.To obtain a MCM node mapping criterion that guarantee an acceptable balance between results accuracy and computational costs, an analysis of the influence of the discretization level obtained by the model was conducted.
Furthermore, to assess the validity of the proposed MCM, three different size IPMSMs, whose characteristics are presented in Table 3, were evaluated.The analysed machines have very different geometries to validate the proposed MCM and establish a node mapping criterion, regardless of the size of the machine.
The load conditions for each machine are presented in Table 4. −60 −79 −34 As Table 5 and Figure 9 illustrate, the three machines operate under heavy load conditions for the considered operating points.In the case of Motor C (Figure 9c), although the stator is not highly saturated, it can be observed that the armature current MMF is comparable to permanent magnet MMF, thereby establishing the magnetic axis almost in the q-axis.Besides, the machine comprises wide slots that behaves as barriers to the magnetic flux and, therefore, increasing the magnetic saturation at the rotor bridges.Altogether, this makes not only the magnetic flux paths more unpredictable but also the cross-coupling effect appreciable.To establish a general node mapping criterion, the presented motor geometries have been evaluated with the aid of the developed MCM and the discretization level presented in Table 6.Each column defines a whole machine node mapping and is defined so that a specific machine's modelled part is studied, i.e., in the case of Stator X, the node mapping parameters that define the stator are modified while the rest of parameters are kept at their lowest possible value.The parametric analysis carried out allows to separately analyse the influence of each model.The Xs in Table 6 refer to values that have been evaluated in each analysis, ranging from 1 to 6.For the parameter controlling the rotor discretization (Rotor X), the Xs values are evaluated in threes due to the complexity of the IPMSM rotor geometry and the difficulties to predict the magnetic flux paths at heavy load conditions.

Parameter StatorX AirgapX BridgeX RotorX
As explained in Section 2.5, the key parameter to determine the machine performance is the magnetic flux linkage, Ψ load,ph .Thus, it is carried out a further analysis on this parameter obtaining.
In Table 7, Ψ load,ph results obtained by FEM are presented.In Table 8, the relative errors between the results of the MCM and FEM are presented for the first and third time harmonics of the Ψ load,ph waveforms.Each row corresponds to results obtained when a MCM node mapping defined in Table 6 is employed with a determined X value.For example, Stator 1 means one value for all the stator parameters, and the other parameters are set to the default value.In Figure 10, the error tendency for the first harmonic is graphically plotted.The discrepancies between the results obtained by the MCM and the FEM are fairly low, especially taking into account the deep saturation of the iron parts for the considered operating points.
To show graphically the great accuracy obtained by the developed MCM, some important magnetic calculated characteristics are presented and compared with obtained by FEM simulation.
In Figure 11, the Ψ load,ph corresponding to Motor A is presented.Its time distribution at an electric cycle for a given stator phase is shown.In Figure 12, the spatial distribution of airgap magnetic flux density along a pole for Motor A at no load operation point is shown.It reflects the accuracy obtained in terms of taking into account the slotting effect, noticing the well calculated magnetic flux density ripple, and also, in terms of the well calculated beginning and ending of the mentioned waveform, which means that the bridge's model works correctly.The MCM node mapping used is the one termed "Stator6", due to the excellent results obtained.Employing MCM magnetic information, the main electromagnetic characteristics at load conditions have been calculated as explained in Section 2 for the three machines.These are presented in Table 9.

Results Discussion and Node Mapping Criterion
When a machine is modelled by a magnetic reluctance network, it is usually not clear which node mapping configuration offers the optimal balance between accuracy and simulation time.With the aim of establishing a node mapping criterion, the results collected in Table 8 and plotted in Figure 10 are examined in greater detail in the following paragraphs:

Stator (Figure 10 a)
As can be observed, the higher is the number of stator nodal elements, the better is the accuracy of the results.The main reason is that, when the number of stator elements increases, so does the number of airgap elements.

Rotor (Figure 10 b)
Although the rotor has a considerable number of nodal elements, the stator is poorly meshed.This means that the linking between airgap and stator is not good enough, especially in terms of taking the slotting effect into account.
For Motor A and Motor B, the number of airgap nodes does not increase as much when N col,PM is increased, unlike the case where the stator is more finely discretized.Nevertheless, for Motor C, the number of airgap nodes is significantly increased, because the number of slots per pole is low, and therefore the effect of stator discretization in the airgap meshing is lower.

Bridge (Figure 10 c)
Because of the high magnetic saturation in the bridge, the number of nodes does not affect the results.

Airgap (Figure 10 d)
The airgap meshing is controlled by the number of existing nodal elements at the stator-airgap border, and the rotor-airgap border.A different number of N row,agap , does not involve a thorough whole machine physical domain modelling.Therefore, the obtained results are invariant.
Additionally, the computational time was measured for the different mapping configurations.In Figure 13, the normalized simulation time is displayed for Motor B using different node densities.For the analysis called BridgeX, the defined MCM has a slightly higher computational cost due to the increasing number of deeply saturated permeance sub-elements.In addition, it should be noted that the high computational cost for the RotorX analysis is due to the fact that the number of rotor elements increases three times faster than parameter X.
The reflected data in Table 10 show the required computation time using FEM and MCM corresponding to Stator1 node mapping.The time data relate to one time step solving process.It should be highlighted that the required computation time for solving MCM shown in Table 10 is up to seven times less than the required time for solving FEM.This together with the fact that, in contrast to FEM, MCM is instantaneously generated, makes our proposed MCM a suitable and fast designing tool.
Based on the accuracy of the results and the required solving time, it can be stated that a good node mapping choice is in the range of configurations Stator3 to Stator6.
In Table 11, the airgap node mapping densities are presented for all three machine configurations under study.The information corresponds to a single airgap row.Various airgap node row details are given: number of nodes per pole pair, average distance between two consecutive nodes, and number of nodes per slot pitch.Next, an analysis to obtain a general node mapping criterion was carried out.

•
The number of nodes per pole pair is increased because of the higher value of X.Additionally, the difference between different geometries is due to the fact that the number of slots per pole are different, and the stator mapping is modelled for each slot-tooth pair.In this regard, setting a number of nodes per pole is not sufficient to establish a general node mapping criterion.

•
The average distance between consecutive nodes depends on the number of airgap nodes and especially on the airgap diameter, i.e., motor size.As can be seen in Motor C, despite having very low values, it does not guarantee excellent results.

•
The number of nodes per slot pitch for Motor A and Motor B with at least a value of 5 is sufficient for obtaining acceptable results.Nevertheless, when analysing rotor mappings, although this ratio is reached, the stator magnetic flux paths are not very well defined.In the case of Motor C, due to the lower number of nodes per pole pair, a value of at least 10 (Stator3) is needed.This value is also reached at rotor analysis, but, as explained previously, the stator node mapping is not very well defined.
The analysis presented here leads to the next node mapping criterion for modelling IPMSM by MCM.
1.The key point is that the number of nodes per slot in the airgap must be at least seven.
As stated in Table 8 and the airgap region nodal information in Table 11, specifically those displayed in the third column (defined with the heading Nodes per slot and pitch), the best results are obtained for MCM node mappings with a ratio value of at least seven nodes per slot in the airgap region.2. There must be at least three equidistant nodes in the stator slot and the tooth span.
To consider the slotting effect and hence obtain accurate results, according to Table 8, the nodal configuration named Stator 3, that guaranties having at least three equidistant nodes in the stator slot and the tooth span, is enough.3.At least three rows of nodal elements show be employed at both rotor yokes.
In consonance with Figure 10, to increase the stator nodal density (Figure 10a), it is fundamental to obtain accurate results, even with the minimum number of rotor nodal rows (N col,PM ), which in this case has been established as three.
By applying these criteria, an adequate connection between the stator and rotor models can be generated, leading to accurate results.

Experimental Validation
With the aim of validating the proposed methodology, a real traction motor with embedded magnets was tested.Its main characteristics are shown in Table 12.The prototype machine has been previously simulated via FEM at rated load.In Figure 14, it can be appreciated that the motor is magnetically saturated, reaching more than 2 T at some points.After applying the node mapping criteria, the main characteristics were obtained and compared with FEM simulations, as shown in Table 13.The operation point is set by defining the same current values for both numerical models.The prototype performance was obtained by standard IEC testing in the test bench shown in Figure 15.In Table 14, the experimental results at rated load conditions are presented.To allow for a fair comparison, the MCM was evaluated at an operating point corresponding to the rated power.Finally, to visually compare the carried out measurements, in Figure 16, the obtained Back electromotive force waveform at tests and MCM is shown.As it can be observed in Table 13, and especially in Table 14 and Figure 16, both experimental and MCM results match largely, showing that the presented approach is suitable for designing IPMSM and predicting its performance even under heavy saturation.

Conclusions
Nowadays, the use of PMSMs is exponentially increasing, with a specific interest in its use at very demanding conditions.Consequently, taking into account the different electromagnetic phenomena that take place inside the machine is crucial to predict in a very precise manner machine performance.
To this aim, in this paper, a general node mapping criterion for modelling highly saturated IPMSM using a magnetic reluctance network has been proposed.A MCM model, based on a magnetic reluctance network, is developed to model V-shaped interior mounted magnet rotors.The proposed MCM model allows selecting the discretization levels for the different machine parts.Furthermore, it not only allows simulating rotor motion but also considers the magnetic cross-coupling effect, the slotting effect and the iron magnetic saturation.
To validate the MCM model, Several V-shaped machines of different sizes and geometries were used together with FEM simulations.The results from these validations were remarkably accurate and efficient, requiring less time to complete the process.Finally, aiming to validate the MCM model with the suggested node mapping criterion, a real IPMSM prototype was tested.The comparison of the obtained results shows a great correspondence, proving the validity of the proposed method to determine highly saturated Interior PMSMs performance.
Overall, the present study expands the field of magnetic circuit models for highly saturated machines.Specifically, it provides the means to evolve within this area towards using the node mapping criterion to apply it to other highly demanded PMSM rotor solutions such as the Spoke type or the multilayer IPMSM.In addition, it would be interesting to analyse the way to optimize the
and embedded in Figure6are shown.

Figure 9 .
Figure 9. Magnetic flux density distribution for the three motors.

Figure 12 .
Figure 12.Airgap magnetic flux density at no load condition (Motor A).

Figure 13 .
Figure 13.Normalized time simulation vs MCM node map.
(a) Calculated by MCM.(b) Measured at test bench.

Table 5 .
FEM measured B max .

Table 6 .
MCM node mapping employed in the analysis.

Table 8 .
MCM node mapping Ψ load,ph relative error values.

Table 10 .
Computational time MCM vs FEM.

Table 13 .
Prototype's performance calculated via proposed MCM and FEM.

Table 14 .
Prototype measured performance at test bench.