Comparative Analysis of the Steady-State Model Including Non-Linear Flux Linkage Surfaces and the Simplified Linearized Model when Applied to a Highly-Saturated Permanent Magnet Synchronous Machine—Evaluation Based on the Example of the BMW i3 Traction Motor

This paper presents a finite element method (FEM)-based model, which describes the magnetic circuit of the BMW i3 traction machine. The model has been reconstructed based on data available in the public domain. The reader is provided with numerical data regarding flux linkage surfaces in dand q-axes, as well as with all the information needed to develop a space-vector model of the machine in steady-state, taking into consideration the non-linearity of the magnetic circuit. Hence, the data of a highly-saturated machine from a renowned product are provided, which can serve as a reference design for research. After that, torque curve and partial load operation points are calculated. Finally, the machine model is linearized and the calculations are repeated with the simplified linearized model. The results from both models are then compared with each other. This comparison is intended to assess the magnitude of the expected inaccuracies, when simplified analytical tools are applied to highly-saturated machines (which are the backbone of automotive electrical drivetrains). It is especially important with regard to preliminary design of electrical drivetrains, as at this stage detailed machine geometry and materials are not known.


Introduction
The electric drive has been gaining popularity in the automotive industry over the past few decades. As automotive applications impose very demanding requirements on drive performance, it has pushed research activities towards the development of highly-efficient, high-power-density, high-torque-density and at the same time cost-effective solutions [1]. As a result, over the years voltage-source inverter fed permanent magnet synchronous machines (PMSM) have become the backbone of not only traction drives, but also drives mounted in auxiliary devices and actuators [1]. In order to meet these strict and often contradicting requirements, PMSM-based drives need to operate in the field-weakening mode and the machines are designed with highly-saturated magnetic circuits [1,2]. Precise analysis of such drives requires complex magnetic field models based on a finite element method (FEM) simulation [3][4][5]. On the other hand, at the preliminary drivetrain design stage, machine dimensions and materials are not known and some simple analytic models are much more desired for this purpose [6].
Most analyses regarding an inverter-fed PMSM operating in the field-weakening region were carried out using analytical formulas basing on a simplified constant parameter model, which assumes linearity of the magnetic circuit [6][7][8][9][10][11]. The main advantage of this approach is that it allows to consider interactions between machine design and power

Materials and Methods
The analysis is performed based on the example of the BMW i3 traction motor. This machine has been chosen for various reasons. First of all, it is an excellent example of a highly-saturated PMSM. Additionally, it is part of a high-end product from a high-tier manufacturer.
The used machine model has been reconstructed based on information available in the public domain. The main source of the data is a tear-down analysis report published by the Oak Ridge National Laboratory (ORNL) [19]. Some additional information about slot geometry, winding configuration and identified materials was provided by prof. Hendershot during his appearance at [20]. System parameters such as DC-bus voltage, maximal current, etc. were given in the presentation from BMW [21].

Machine Geometry and Winding Configuration
The machine model was created in the FEM electromagnetic field simulation software ANSYS Maxwell 2D (2018.1.0) [30]. The modeled geometry along with the most important dimensions are presented in Figure 1. These dimensions were published in the ORNL report [19]. As exact rotor geometry information was not published, it was reconstructed based on a photograph from the same report. Slot geometry is shown in Figure 2. The main dimensions in this picture were presented in [20]. The remaining dimensions were reconstructed based on a picture included in [20].
The winding configuration of the machine was presented in [20] as well. It consists of six Y connected three-phase windings with isolated neutral points, which are connected in parallel. Therefore, there is one Y system per pole pair. A single phase winding of each system is placed in four slots, with nine turns per coil. The exact placement of coil sides in particular slots is presented in Figure 3. A, B and C are the phase names. The symbol '+' means the forward direction coil side and the symbol '-' means the return direction coil side.

Definition of the Used Materials
The materials used in the machine were identified by ORNL and provided in [20]. From the magnetic circuit simulation point of view, the most important are parameters of the used steel and magnets.
It has been identified that the rotor core, stator yoke and stator teeth are all made of the same material, i.e., M-19 29 Ga steel. The magnetization curve of this steel shown in Figure 4 has been reconstructed based on the characteristics given in [31].
The rotor magnet material has been identified to be Neodymium Iron Boron 38/23. According to [32], the characteristics of this material are: magnetic coercivity H c = 955 kA/m; magnetic remanence B r = 1.24 T.

Model Summary
For the readers' convenience, all important machine parameters have been summarized in the form of tables in Appendix A. Mechanical dimensions are provided in Table A1, winding information in Table A2 and material information in Table A3.
Please note that the presented FEM model has also been included in the Supplementary Materials attached to this paper (in the Ansys Maxwell file format).
For readers using other FEM simulation software, the machine geometry has been additionally exported in the form of a .dxf file. This geometry together with the information from previous subsections is sufficient to fully recreate the model. Figure 5 shows exemplary flux density simulation results obtained with the model. It is important to note that the used steel saturates rapidly above approximately 1.25 T (see Figure 4). The first field plot (Figure 5a) shows the state when no currents are fed into the windings. In this state, the magnetic circuit is saturated at magnet sides. This is a necessary feature of every interior permanent magnet synchronous machine (IPMSM) design [1], as it provides very high reluctance to the leakage flux path around magnet sides and therefore forces the magnetic flux out of the rotor through the air gap. It is worth mentioning that even in this current-less state the rotor yoke is saturated at a flux return path (which indicates relatively high core material utilization). The second plot (Figure 5b) shows the state when currents resulting in the peak torque of the drive are fed into the windings (it is explained later in this paper how the values of these currents were obtained). It can be seen that these currents cause significant saturation in the stator yoke, the stator teeth, as well as bulk rotor saturation.

FEM Simulation Results
The FEM simulation was carried out for many different combinations of the dand q-axis currents and for one full electrical period. It should be taken into consideration that the rotor of this machine is segmented in six parts in the axial direction (see the report [19]). These segments are rotated consecutively by one mechanical degree providing rotor skewing, which helps to partly eliminate adverse spatial harmonics in torque and induced voltage. It was considered in the analysis in such a way that the simulation was carried out for six different initial rotor positions and the results were averaged.
Please note that all the results (after averaging due to rotor skewing) are provided to the reader in Supplementary Materials attached to this paper. The attached spreadsheet includes the following signals as functions of mechanical angle: machine torque, d-axis flux linkage, q-axis flux linkage and induced voltages in all three phases.  Tables A4 and A5, respectively. Torque values can be found in Table A6. The peak to peak torque ripple values were calculated as well (see Table A7). Values for intermediate current operation points were obtained via cubic interpolation and plotted in the form of three dimensional surfaces in Figure 6. Severe magnetic circuit saturation due to an increase of the q-axis current component can be clearly seen in the q-axis flux linkage surface ( Figure 6b). On the other hand, the d-axis flux linkage surface is relatively linear (Figure 6a), although it bends slightly downwards with the rising q-axis flux linkage value in a low d-axis current range. This is due to the so called cross-saturation effect caused by the fact that both dand q-axis fluxes share some common path in the magnetic circuit [18]. This effect is also visible in the q-axis flux linkage surface, as the d-axis de-magnetization (i.e., rise in the de-magnetizing d-axis current magnitude) causes a slight increase in the q-axis flux linkage value. All the above findings indicate that the modeled machine exhibits significant iron core saturation in its feasible operation range. Hence, it has proven to be a very good candidate for comparison with the linearized model, since the obtained discrepancy between the results calculated with these two models can be expected to be close to its realistic extreme values.
The model was validated based on the locked rotor test results published in the ORNL report [19] (see Figure 7). During such a test the rotor is locked and currents of constant amplitude and various space-vector angles (with respect to d-axis) are fed into the windings. Importantly, the machine torque contains spatial harmonics and thus the value measured in the locked rotor test depends on the position in which the rotor has been locked. In order to take it into consideration, the simulation values used for the comparison in Figure 7 are the values averaged over one electrical period and the expected torque ripple is shown in the form of shaded areas. For this validation the torque and the torque ripple values presented in Figure 6c,d were used. The values corresponding to the measurement points were calculated using cubic interpolation. It should be considered that the rotor geometry of the machine is relatively complex and it was reconstructed based only on a photograph. Hence, it is obvious that some inaccuracies of the model should be expected. A difference between torque values calculated with the FEM model and the measurement results is especially visible at angles around 90 degrees (i.e., for low d-axis current). In the authors' opinion, the obtained accuracy is sufficient to positively validate the presented model.

Calculations Based on the Non-Linear Model
The torque vs. speed curve can now be calculated using the following non-linear steady-state model [18]: where u d and u q are dand q-axis terminal voltage space-vector components (V); i d and i q are dand q-axis current space-vector components (A); R is stator phase resistance (Ω); ω el is electrical rotor speed (rad/s); ψ d (i d , i q ) and ψ q (i d , i q ) are dand q-axis flux linkages (Wb), which are non-linear functions of the current space-vector components (i.e., i d and i q ); p is the number of pole pairs; T M is machine torque (Nm). The torque maximizing operation points were calculated based on (1) and the values listed in Table 1 using a numerical search algorithm created in MATLAB (R2017b) software [33]. The results are shown in Figure 8a, where the filled dots represent the obtained operation points for the consecutively rising speed. The same points can be identified on the torque vs. speed curve shown in Figure 8b. The first one corresponds to the peak torque in the base operation speed range and lies in a place where the possibly uppermost torque iso-line (blue lines in Figure 8a) is tangential to the current limitation circle (red circle in Figure 8a). The flux density plot in Figure 5b was obtained for exactly this operation point, which can be described with the following current space-vector component values: It is common knowledge that there exist two general drive classes regarding their behavior during field-weakening operation. The first is the finite maximal speed drive class. In such a case, there is only one field-weakening operating mode, i.e., current and voltage limited operating mode. During such operation, the current space-vector locus follows a current limitation circle up to the maximal speed of the drive. The second class is the infinite maximal speed drive class. In this case, there exists an additional operating mode above some threshold speed. This is a voltage limited operation according to a so-called maximal torque per voltage (MTPV) control strategy. During such an operation, the current space-vector locus follows a path from the current limitation circle into the so-called machine characteristic current point. All these characteristics are explained in detail in [6,10].
It can be identified in Figure 8a that the analyzed drive exhibits only one fieldweakening operation mode, i.e., current and voltage limited operation along the current limitation circle path. This operation mode is maintained up to the maximal speed of this drive (i.e., 11,400 rpm) and the MTPV operation region has not been identified (which means that this is the finite maximal speed drive according to [6,10]).

Calculations Based on the Linearized Model
The goal of this paper was to quantify the influence of iron core saturation and voltage drop across the winding resistance on the results obtained with the following simplified model (for the details of this model please refer to [6,10]): where L d_lin and L q_lin are dand q-axis inductances calculated for some particular linearization point (H); ψ PM_lin is permanent magnet flux linkage calculated for some particular linearization point (Wb). In order to do so, the torque vs. speed curve for the BMW i3 traction drive should be now calculated with analytical formulas derived using this model, which can be found in [6]. The results can be then compared with these obtained using the previously introduced more complicated model (1). For convenience, in the course of this paper the following names are going to be used, in order to distinguish between the two models: • 'Non-linear model'-the space-vector model in steady-state defined with Equation (1), which includes voltage drop across stator resistance and non-linear flux linkage surfaces. • 'Linearized model'-the simplified model defined with Equation (3), which neglects voltage drop across stator resistance and assumes constant machine parameters (hence the term 'linearized').
As the linearized model uses constant machine parameter values, the non-linear flux linkage surfaces obtained in FEM simulation (see Figure 6) should be linearized at some operation point. In the presented analysis, the peak torque operation point (2) was chosen. The linearized machine parameters can be calculated from the flux linkages definition: as follows: The obtained linearized system parameters have been summarized in Table 2. These data follow the nomenclature introduced in [6]: where ε is the machine saliency factor; k ch is the drive characteristic factor; i max is the maximal phase current of the drive (A); i ch is a machine characteristic current, i.e., d-axis current magnitude needed to reduce the d-axis flux linkage to zero (A). It is a well-known fact that the field-weakening performance of an inverter-fed PMSM drive depends on the relationship between the motor characteristic current i ch and the maximal current of the drive i max [8,10,13]. In [6], the authors proposed to describe this relationship with a single variable, i.e., the drive characteristic factor k ch . It was also derived that the value of this factor can deliver information about the necessity to over-size the drive inverter power rating.
For k ch ≤ 1, the drive has a finite maximal speed and the ratio of this speed to the base speed rises with the value of this factor. The peak power value of drives from this class is equal to the volt-ampere rating of the inverter needed to operate these drives within their specifications. It means that for this class the power rating of the power electronic converters does not need to be over-sized.
For k ch > 1, the drive has an infinite maximal speed and the true constant power speed region rises when this factor is increased. Unfortunately, together with the value of k ch the ratio between the volt-ampere rating of the power electronics inverter needed to operate the drive and the peak power at the machine shaft also rises. Hence, in this case the inverter power rating needs to be over-sized. A detailed discussion regarding this topic can be found in [6].
For the reasons explained above, it is common practice to design the drives in such a way that the characteristic factor has a value possibly close to the unity. It provides a good trade-off between field-weakening performance and power converter sizing. The analyzed BMW i3 drive is a perfect example of this design trend. Please note that the obtained drive characteristic factor value for this drive (see Table 2) is slightly smaller than the unity, but very close to it. As the drive belongs to the finite maximum speed class (i.e., k ch < 1), there should be no MTPV operation region, which matches the results obtained with the non-linear model (see Figure 8). Now, the torque vs. speed curve can be calculated based on the linearized model using analytical formulas. Please note, that different equivalent forms of these equations can be found in the literature [6,8,10]. The solution derived in [6] is going to be rewritten here for the reader's convenience. It bases on the following normalized quantities in a per-unit system: where the superscript (• * ) denotes the per-unit quantity (p.u.); i is current (A); ω el is electrical angular speed (rad/s); T is torque (Nm); u max is the maximal phase voltage (V); ω N is the base electrical angular speed for the normalization (rad/s); T N is the base torque for the normalization (Nm). The following formulas allow to calculate per-unit values of currents, torque and speed, which can be converted into physical quantities with (7). The per-unit peak torque in the base speed region can be calculated with: The per-unit transition speed between the base speed region and the field-weakening speed region equals: The per-unit torque in the field-weakening speed operation range can be expressed as a function of the per-unit speed as:

Comparison of the Results-Maximal Torque
The results of the maximal torque vs. speed calculations for both methods are shown in Figure 9 (all data for the non-linear model are marked with black and for the linearized model with gray). It can be seen that both torque surfaces diverge significantly from each other in the high q-axis current region (see Figure 9a). On the other hand, matching between the surfaces is relatively good in the low q-axis current region. It can be observed, that operation points in the field-weakening speed range (i.e., above ca. 4000 rpm) lie in this quasi linear region, hence the torque vs. speed curves obtained with the two methods match very well in this speed range (see Figure 9b). On the other hand, a vast part of the current limitation boundaries (solid quasi quarter-circular curves on the surfaces) lies in the saturated region of the torque surface (see Figure 9a). It is a well-known fact [6][7][8][9][10][11][12] that the peak torque operation point lies somewhere on this boundary and since both surfaces diverge significantly in this region, completely different peak torque operation points were identified with both methods. As a result, there is a relatively big difference in the peak torque value obtained with both models: 279.7 Nm with the linearized model vs. 258.2 Nm with the non-linear model.

Comparison of the Results-Partial Load Operation
Based on the above-mentioned findings, it is then interesting to obtain linearized model inaccuracy in calculating the torque value for partial load operation. For this reason, the torque vs. speed curves were calculated with both models once more, but this time for many different values of the maximal current. The results can be found in Figure 10b. It can be noticed that the differences between torque values calculated with both models in the base speed region are much bigger than in the field-weakening speed region. Hence, a more detailed analysis of this operational region is needed. The calculation results for this case are presented in Figure 11. The current space-vector loci for the peak torque calculated with different current limitation values are presented in Figure 11a. Again, the solid quasi quarter-circular curves on the surfaces indicate current limitation curves for different current values. A corresponding torque vs. phase current magnitude plot is shown in Figure 11b. The curves were calculated for values up to the 650 A (which is a 115% of the feasible operation range of this drive) for illustrative purposes. The blue filled dots indicate the actual maximal phase current limitation of the drive (i.e., 565.7 A).

Discussion
It can be identified that torque values in the field-weakening speed region match relatively well for both models. The current space-vector loci calculated with both methods follow exactly the same path in the dq-plane (see Figure 10a). These paths are current limitation circles for the consecutive maximal current values. This is going to be the case for every finite speed drive, as for this drive class, this path is predefined. It should be mentioned that these results can not be extended to the infinite speed drive class, because the current locus path for these drives during the MTPV operation is not predefined. It can be then expected that the results obtained with both models can diverge from each other much more, as in the case presented here.
There is also one more important observation, which should be pointed out clearly. Even though the torque vs. phase current magnitude curves for each model match relatively well with others in the partial load operation range, the corresponding current operation points in the dq-current plane are very different. It means that the linearized model is incapable of calculating proper current space-vector components for the maximal torque per current (MTPC) operation (see [6,10]).
In the authors' opinion, the results obtained with the linearized model match surprisingly well these obtained with the non-linear one (it should be remembered that the analyzed example is an extreme case of a highly-saturated machine). These results are rather unexpected, as it seems to be common belief that the linearized model is not feasible for calculating the torque performance of highly-saturated machines. Even though peak torque for the actual maximal current was calculated with a 21.5 Nm (8%) error, torque values for the partial load and field-weakening operation are still relatively accurate.
Having some basic impression of the extent to which the magnetic circuit of the designed machine is going to be utilized (in terms of material magnetic saturation), some error margin can be assumed and the linearized model can be still successfully applied for the torque vs. speed curve estimation. It is important to emphasize that the conclusions drawn in this paper should be limited to the finite maximal speed drive class. Only for this class the current locus path during the field-weakening operation is predefined to go along the current limitation circle and this is probably the cause of such good performance of the linearized model in the field-weakening speed region.

Conclusions
Complex information about the BMW i3 traction machine has been provided, i.e., mechanical dimensions, winding layout information and material definition (see Tables A1-A3). The presented FEM model and its geometry are included in the Supplementary Materials. The FEM simulation results including spatial harmonics information are given too.
Torque and flux linkage values averaged over one electrical period are provided as well (see Tables A4-A6). These numerical values, together with data from Table 1 fully define the space-vector model of this drive during steady-state operation, which also takes into consideration non-linearity of the magnetic circuit. Along with the references [19][20][21], the reader is provided comprehensive information about a high-end industrial product, which can be used as a reference design during research studies.
It has been shown that the presented motor is an extreme example of a highly-saturated IPMSM. Based on the presented results, it can be concluded that to some extent the simplified linearized machine model can be successfully applied for the torque vs. speed curve estimation of a highly-saturated machine. Nevertheless, some proper error margin should be taken into consideration.
It seems obvious that using exact information about the magnetic circuit non-linearity is the correct way to obtain accurate results. Nevertheless, in the authors' opinion, it was very valuable (from both practical point of view and out of pure scientific curiosity) to explore the magnitude of the inaccuracies obtained, when the linearized model is applied instead.

Data Availability Statement:
The data presented in this study are available partly in article and partly in Supplementary Material.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: