High Accuracy Modeling of Permanent Magnet Synchronous Motors Using Finite Element Analysis

: Permanent magnet synchronous machines (PMSMs) have garnered increasing interest because of their advantages such as high efﬁciency, high power density, wide speed range, and fast dynamics. They have been employed recently in several industrial applications including robotics and electric vehicles (EVs). However, PMSMs have highly nonlinear magnetic characteristics, especially interior PMSMs, due to the existence of reluctance torque. Nonlinearity complicates not only machine modeling but also control algorithms. An accurate machine model is the key aspect for the prediction of machine performance as well as the development of a high-performance control algorithm. Hence, this paper presents an accurate modelling method for PMSMs. The proposed model method is applicable for all PMSMs, even multiphase machines. This paper considers a fractional slot concentrated winding 12/10 interior PMSM (IPMSM) for this study to demonstrate the effect of magnetic saturation and special harmonics. The developed model considers accurately the magnetic saturation, mutual coupling, spatial harmonics, and iron loss effects. It utilizes ﬁnite element analysis (FEA) to estimate the precise magnetic characteristics of IPMSM. The ﬁnite element model is calibrated precisely using experimental measurements. The iron losses are estimated within the simulation model as d- and q-axes current components. The model accuracy is validated experimentally based on a 12/10 IPMSM prototype.


Introduction
Because of their superiorities permanent magnet synchronous machines (PMSMs) have garnered increasing interest. They have high efficiency, high power density, a wide range of operating speeds, small size and weight, and low noise [1][2][3]. Recently, interior PMSMs (IPMSMs) showed the best overall performance for electric vehicles (EVs) [4,5]. IPMSMs can easily fulfil vehicle requirements with proper torque control [6,7]. However, IPMSMs are well known for their nonlinearity and torque ripple due to their magnetic saturation, spatial harmonics, and reluctance torque [8].
Developing and optimizing different control strategies of IPMSMs for high performance drives rely basically on modelling accuracy. The conventional machine models (CMMs) were introduced using fixed machine parameters regardless of the effects of saturation, spatial harmonics, and iron loss. These models cannot represent accurately the behavior of IPMSMs [3]. They can be used for initial machine designs. Therefore, several research studies have been done involving the production of high-fidelity machine models. The circuit-field coupled co-simulation is presented in [9][10][11]. This technique is very time-consuming due to the involvement of numerical finite element (FE). In [12], the relationship between flux linkage and currents are fitted to be used within the simulation in order to avoid the FE computation in simulation. However, the fitting is a very complicated and time-consuming process. In [13], the modelling in the abc-reference frame is introduced using N-dimensional look-up tables (LUTs). The abc-model allows one to include both geometry and saturation-induced space harmonics naturally. The abc-model can provide relatively high accuracies, while the complex expressions are circumvented. In [14], conventional shallow neural networks (CSNN) are applied to fit torque-current characteristics. Due to the high nonlinearity of torque-current characteristics as well as the limited capability of CSNN to fit complex nonlinear systems, the torque and the significant harmonics cannot be accurately modelled. To consider the torque nonlinearity and the spatial harmonics, deep neural network are developed in [15]. However, the training of intelligent techniques requires high skills and a large number of given data. Moreover, measured samples are required to train the network or generate the rules.
Considerable interest has been given to the analytical models for IPMSMs. Modelling and quantifying the source of torque ripple as part of the dq-model are proposed in [16]. In [17,18], accurate estimation of average torque, including nonlinearity effects, is presented. In [8], with the consideration of magnetic saturation, spatial harmonics, cross-coupling, and temperature effects, the torque and the voltage equations are derived. However, these analytical models are very complicated for real-time implementation. In [19], an indirect analytical method is presented for V-shaped IPMSM. It transfers the rotor into an equivalent linear one in a polar coordinate system to analyze the magnetic field and predict cogging torque. However, a considerable error is obtained in predicted air-gap flux density. Moreover, the torque calculation error gradual increases with rise of current. In [20], the flux-linkage map of IPMSM is identified by experimental tests in motoring action. However, although this method can consider the magnetic saturation and spatial harmonics, it has a very complicated and time-consuming procedure that involves several measurements.
For accurate modelling, the iron losses must be included as they affect output torque. In [21,22], the iron loss analysis was derived based on FEA. In [23], the numerical calculation of iron loss' resistances based on FEA are presented. However, the modelling of iron loss would be very time-consuming using FEA. Hence, analytical methods are used. In [24], the calculation of iron loss of PMSM's stator considering DC-biased magnetic induction is introduced. Moreover, in [25,26], the iron losses are expressed in the form of d-and q-axes flux linkages or currents. However, the derivation of iron loss' coefficients still remained unsolved. Moreover, the effect of iron loss on electromagnetic behavior of IPMSMs have to be considered and analyzed.
This paper develops a high-fidelity and accuracy model for PMSMs. The proposed model fits all types of PMSMs. It represents accurately the performance of real machines by precise consideration of magnetic saturation, mutual coupling, spatial harmonics, and the iron loss effects. First, the FEA is employed to generate the magnetic characteristics of the tested machine. Then, these characteristics are validated experimentally. The experimental measurements aim basically to have a precise and calibrated finite element (FE) model that represents a real machine. The experimental measurements are very simple, and they include the measurement of PM flux linkage and the calibration of the BH curve of the magnetic steel material. Once, the FE model is calibrated precisely, it is used to generate full data of the machine. Second, the FE data are processed and rearranged to produce a simple and computationally efficient MATLAB Simulink model. Third, the effect of iron losses is achieved analytically and integrated within the MATLAB Simulink model. This paper considers a fractional slot concentrated winding 12/10 IPMSM. The experimental verifications are included to prove model accuracy and flexibility.

Conventional Model of IPMSM
The dynamic model of IPMSM in a d-q rotor reference frame is as follows [27][28][29]: where v d and v q , i d and i q , λ d and λ q , and L d and L q are dand q-axes components of voltage, current, flux-linkage, and inductance, respectively. λ PM depicts the flux-linkage of permanent magnets. R s is the stator resistance. ω e is the electrical angular speed. The electromagnetic torque (T e ) of IPMSM has two components, namely, the reluctance torque and the magnet torque, as described in Equation (3).
where p is the pole pairs. The mechanical equation is described by Equation (4).
where J is the inertia, ω m is the angular speed, B is the frictional coefficient, and T L is the load torque. Figure 1 shows the implementation of the conventional machine model (CMM) using Equations (1)-(4). The CMM uses a constant machine parameter. Hence, it neglects the magnetic saturation, mutual coupling, spatial harmonics, and iron loss effect. The CMM considers only the fundamental components and neglects harmonic fields that are the result of magnetic saturation, slotting, and variation of permeance with rotor position. The inaccurate model affects performance analysis and accurate prediction of drive system performance. Moreover, it compromises control quality in maximum torque per ampere (MTPA) and field weakening (FW) operations. Hence, an accurate model that includes magnetic saturation, mutual coupling, spatial harmonics, and iron loss effects is necessary.
where vd and vq, id and iq, λd and λq, and Ld and Lq are d-and q-axes components of voltage, current, flux-linkage, and inductance, respectively. λPM depicts the flux-linkage of permanent magnets. Rs is the stator resistance. ωe is the electrical angular speed. The electromagnetic torque (Te) of IPMSM has two components, namely, the reluctance torque and the magnet torque, as described in Equation (3).
where p is the pole pairs. The mechanical equation is described by Equation (4).
where J is the inertia, ωm is the angular speed, B is the frictional coefficient, and TL is the load torque. Figure 1 shows the implementation of the conventional machine model (CMM) using Equations (1)-(4). The CMM uses a constant machine parameter. Hence, it neglects the magnetic saturation, mutual coupling, spatial harmonics, and iron loss effect. The CMM considers only the fundamental components and neglects harmonic fields that are the result of magnetic saturation, slotting, and variation of permeance with rotor position. The inaccurate model affects performance analysis and accurate prediction of drive system performance. Moreover, it compromises control quality in maximum torque per ampere (MTPA) and field weakening (FW) operations. Hence, an accurate model that includes magnetic saturation, mutual coupling, spatial harmonics, and iron loss effects is necessary.
To consider only the magnetic saturation, the d-and q-axes inductances (Ld, Lq) are presented as functions of id and iq. The obtaining of such relationships is not an easy task; several measurements must be conducted, which complicate the idea of obtaining a high-accuracy model for an IPMSM [30].

The Proposed Machine Model
For real IPMSMs, the fluxes (λd and λq) and the electromagnetic torque (Te) are functions of currents (id and iq) and rotor position (θi). These relationships can be written To consider only the magnetic saturation, the dand q-axes inductances (L d , L q ) are presented as functions of i d and i q . The obtaining of such relationships is not an easy task; several measurements must be conducted, which complicate the idea of obtaining a high-accuracy model for an IPMSM [30].

The Proposed Machine Model
For real IPMSMs, the fluxes (λ d and λ q ) and the electromagnetic torque (T e ) are functions of currents (i d and i q ) and rotor position (θ i ). These relationships can be written as given by Equation (5). Such an equation considers the saturation and spatial harmonics. Moreover, the torque components as well as the cogging torque are also considered. The proposed machine model is developed as shown in Figure 2 considering the given fact of Equation (5). The inputs are the voltages v d , and v q . The flux-linkages (λ d and λ q ) are estimated using Equation (6). The i da and i qa currents are obtained from the flux inverse model that is given by Equation (7). The motor torque is obtained as a function of currents (i da , i qa ) and rotor position θ i . The iron losses are represented by the currents i d-fe and i q-fe . The estimation of iron losses within the model is given in a later section. Finally, the equation of mechanical motion is used to estimate the rotor position and speed.
The degree of modeling accuracy depends on calculation accuracy of ma characteristics (λd (id, iq, θi), λq (id, iq, θi), and Te (id, iq, θi)) for the tested machine. IPMSM possess high nonlinear magnetic characteristics, the estimation of λd (id, iq (id, iq, θi), and Te (id, iq, θi) is not straightforward. Therefore, the FEA is emplo estimate the magnetic characteristics for the tested 12/10 IPMSM prototype. Figure 3 shows the machine structure. It is a 12 slot, 10 magnet pole, 3 IPMSM. The rotor magnets are the V-shaped type. The stator has double layer wi for better torque ripple reduction. The nominal parameters and dimensions of I are shown in Table 1. The degree of modeling accuracy depends on calculation accuracy of magnetic characteristics (λ d (i d , i q , θ i ), λ q (i d , i q , θ i ), and T e (i d , i q , θ i )) for the tested machine. As the IPMSM possess high nonlinear magnetic characteristics, the estimation of λ d (i d , i q , θ i ), λ q (i d , i q , θ i ), and T e (i d , i q , θ i ) is not straightforward. Therefore, the FEA is employed to estimate the magnetic characteristics for the tested 12/10 IPMSM prototype. Figure 3 shows the machine structure. It is a 12 slot, 10 magnet pole, 3 phase IPMSM. The rotor magnets are the V-shaped type. The stator has double layer windings for better torque ripple reduction. The nominal parameters and dimensions of IPMSM are shown in Table 1.

considered.
The proposed machine model is developed as shown in Figure 2 considering the given fact of Equation (5). The inputs are the voltages vd, and vq. The flux-linkages (λd and λq) are estimated using Equation (6). The ida and iqa currents are obtained from the flux inverse model that is given by Equation (7). The motor torque is obtained as a function of currents (ida, iqa) and rotor position θi. The iron losses are represented by the currents id-fe and iq-fe. The estimation of iron losses within the model is given in a later section. Finally, the equation of mechanical motion is used to estimate the rotor position and speed.

=
, , ; = , , ; = , , The degree of modeling accuracy depends on calculation accuracy of magnetic characteristics (λd (id, iq, θi), λq (id, iq, θi), and Te (id, iq, θi)) for the tested machine. As the IPMSM possess high nonlinear magnetic characteristics, the estimation of λd (id, iq, θi), λq (id, iq, θi), and Te (id, iq, θi) is not straightforward. Therefore, the FEA is employed to estimate the magnetic characteristics for the tested 12/10 IPMSM prototype. Figure 3 shows the machine structure. It is a 12 slot, 10 magnet pole, 3 phase IPMSM. The rotor magnets are the V-shaped type. The stator has double layer windings for better torque ripple reduction. The nominal parameters and dimensions of IPMSM are shown in Table 1.

Finite Element Analysis (FEA) of IPMSM
Despite the high accuracy of FEA compared to analytical methods, it needs the exact material properties and geometrical dimensions of the tested machine [31]. The main obstacle appears with the material properties, especially the PMs and magnetic steel properties (BH curve). If the material properties are given precisely, the FEA can provide very high accuracy of magnetic characteristics.
To ensure the inputted data to the FEA model represent accurately the real machine, this paper introduces a calibration method for the FEA model. The calibration is achieved based on the experimental measurements for the tested 12/10 IPMSM prototype. The experimental calibration is very important, as it can consider the manufacturing effects on both magnets and magnetic steel materials. It also can include the manufacturing tolerances and the physical changes that are included due to the manufacturing process.

Calibration of FEA Model
The aims of calibration is to obtain the real material properties for the tested IPMSM. These material properties will be inputted to the FEA model. Any error in material properties will lead to a calculation error in the magnetic characteristics of IPMSM. Therefore, the calculation accuracy of the FEA model depends basically on the accuracy of its inputted data.
The properties of permanent magnets (PMs) rarely coincide with the available data in manufacturers' datasheets. This is basically due to the manufacturing quality and tolerances, especially for low-cost motors. Probably, it will be decreased slightly less than the manufacturer datasheets [32]. Hence, the actual PM flux linkage must be measured experimentally to represent the real machine precisely.
The PM flux linkage (λ PM ) can be measured experimentally when the machine is working as a generator under no load. The no-load back EMF voltages are measured and recorded, preferably at high speeds. The higher speed is preferred to have a higher voltage on stator windings. The measurement errors in the case of a higher voltage range are lower compared to the case of lower voltage. Hence, high speed is preferred to reduce the measurement error in phase voltages and hence to estimate accurately the flux linkage of permanent magnets. λ PM can be estimated using Equation (8). Due to the concentrated windings configuration, the fundamental component of phase voltage must be considered. Hence, Fast Fourier Transform (FFT) is conducted for the phase voltage, as seen in Figure 4.
First, the machine is operated as a generator with the help of a dc motor at a speed of 1708 r/min. The high speed (greater than rated) helps to reduce voltage error and hence to improve the measurement accuracy. The window of FFT is chosen as 5 cycles to ensure accurate analysis, as seen in Figure 4. The analysis is conducted for three phases (A, B, and C) to extract the average value of measured λ PM that represents the real machine properly. The differences appear to be very small and can be ignored, as illustrated in Table 2.
where f is the frequency of back EMF voltage.
Mathematics 2022, 10, 3880 6 of 20 machine properly. The differences appear to be very small and can be ignored, as illustrated in Table 2.
where f is the frequency of back EMF voltage.  Lately, due to the huge interest of EVs, several magnetic steel materials have been introduced. These materials have gone through several tests, including the measurement of their BH curves. The tested motor is made from 35WW270 magnetic steel material. Its measured BH curve is available online. This curve fits accurately with the real machine. Hence, no calibration is needed for the BH curve of magnetic steel.
In case the BH curve data are not available in an accurate enough manner, the machine is operated as a motor with known currents id and iq. Then, the developed motor torque is measured and recorded. After that, the measured torque is compared to the estimated torque from FEA under the same current values (id and iq). The difference, if it exists, is adjusted by the BH curve of magnetic steel.
For now, the accurate coercivity of PM material is obtained. Moreover, the measured BH curve is also available online. Now the FEA model is ready to provide the precise magnetic characteristics for the tested machine. A 2D FE model for the tested IPMSM is illustrated in Figure 5. The motor dimensions are given in Table 1. The model is built in finite element method magnetics (FEMM). The end effects of coils can be ignored due to the considerable axial length of the machine; hence, 3D analysis, which is very time consuming, is not necessary [33]. For the 2D FE model, the recent programs for FEA make the calculation much easier and speed up computations by static magnetic field analysis. The following assumptions are made: a zero magnetic vector potential for the outer surface of the stator; a constant airgap length; isotropic rotor and stator materials; concentric stator and rotor poles;  The measured average PM flux linkage (λ PM ) is 135.95 mWb. The FEA model is then calibrated to produce the same flux linkage (λ PM = 135.95 mWb). Hence, based on this measured value of λ PM , the coercivity of PM in the FEA model is set to 848,000 A/m instead of 969,969 A/m from the manufacturers' datasheets. The reduction ratio is 12.23%.
Lately, due to the huge interest of EVs, several magnetic steel materials have been introduced. These materials have gone through several tests, including the measurement of their BH curves. The tested motor is made from 35WW270 magnetic steel material. Its measured BH curve is available online. This curve fits accurately with the real machine. Hence, no calibration is needed for the BH curve of magnetic steel.
In case the BH curve data are not available in an accurate enough manner, the machine is operated as a motor with known currents i d and i q . Then, the developed motor torque is measured and recorded. After that, the measured torque is compared to the estimated torque from FEA under the same current values (i d and i q ). The difference, if it exists, is adjusted by the BH curve of magnetic steel.
For now, the accurate coercivity of PM material is obtained. Moreover, the measured BH curve is also available online. Now the FEA model is ready to provide the precise magnetic characteristics for the tested machine. A 2D FE model for the tested IPMSM is illustrated in Figure 5. The motor dimensions are given in Table 1. The model is built in finite element method magnetics (FEMM). The end effects of coils can be ignored due to the considerable axial length of the machine; hence, 3D analysis, which is very time consuming, is not necessary [33]. For the 2D FE model, the recent programs for FEA make the calculation much easier and speed up computations by static magnetic field analysis. The following assumptions are made: a zero magnetic vector potential for the outer surface of the stator; a constant air-gap length; isotropic rotor and stator materials; concentric stator and rotor poles; neglected core losses, MMF space harmonics, end-winding, and skin effects; and z-directed current density and magnetic vector potentials. Special attention is directed to the air-gap flux calculation, and a smaller mesh size is used for it. The FEA has 24,838 nodes and 49,292 elements. Figure 5 shows the distribution of flux density as well as the flux lines that are estimated from the FEA model. As noted, the highest flux densities are seen in rotor rips and stator pole shoes. These parts are in deep saturation. neglected core losses, MMF space harmonics, end-winding, and skin effects; an directed current density and magnetic vector potentials. Special attention is direct the air-gap flux calculation, and a smaller mesh size is used for it. The FEA has 2 nodes and 49,292 elements. The results of FEA are obtained as functions of currents id and iq as well as θi currents id and iq are changed from −25 to 0 A (26 points) and from 0 to 25 A (26 po respectively. The current step is 1A. θi is changed from 0° to 72° insteps of 1 mech points). Therefore, it requires (26 * 26 * 73 = 49,348) solutions to generate the magnetic characteristics for the tested machine motor. The symmetry of the mag machine structure helps greatly to reduce the time and/or number of solutions eith plotting part of the machine or by choosing part of the rotor rotation.

Results of FEA
The inductances Ld and Lq are calculated using Equations (9) and (10), respectiv Figure 6 shows the flux-linkages (λd (id, iq, θi), λq (id, iq, θi)), torque Te (id, iq, θi) inductances (Ld (id, iq, θi), Lq (id, iq, θi)) as functions of id, iq at θi = 0. As seen, the linkages, torque, and inductances show highly nonlinear relations with currents. As in Figure 6, the d-and q-axes inductances (Ld, Lq) decrease with the increase of q current iq over the full range of d-axis current id. The results of FEA are obtained as functions of currents i d and i q as well as θ i . The currents i d and i q are changed from −25 to 0 A (26 points) and from 0 to 25 A (26 points), respectively. The current step is 1A. θ i is changed from 0 • to 72 • insteps of 1 mech • (73 points). Therefore, it requires (26 * 26 * 73 = 49,348) solutions to generate the full magnetic characteristics for the tested machine motor. The symmetry of the magnetic machine structure helps greatly to reduce the time and/or number of solutions either by plotting part of the machine or by choosing part of the rotor rotation.
The inductances L d and L q are calculated using Equations (9) and (10), respectively. Figure 6 shows the flux-linkages (λ d (i d , i q , θ i ), λ q (i d , i q , θ i )), torque T e (i d , i q , θ i ), and inductances (L d (i d , i q , θ i ), L q (i d , i q , θ i )) as functions of i d , i q at θ i = 0. As seen, the flux linkages, torque, and inductances show highly nonlinear relations with currents. As seen in Figure 6, the dand q-axes inductances (L d , L q ) decrease with the increase of q-axis current i q over the full range of d-axis current i d .

Production of MATLAB Simulink Model
The proposed machine model that is given in Figure 2 consists of three main parts: first, the torque estimation T e (i d , i q , θ i ), which is a straightforward task, as the obtained data from FEA can be stored in the form of a 3D LUT; second, the representation of iron losses by the currents i d-fe and i q-fe , which is included in a later section; and third, the estimation of currents i d (λ d , λ q , θ i ) and i q (λ d , λ q , θ i ). As noted, this is an inverse solution for the obtained data from FEA. The FEA gives the fluxes as functions of currents and positions. The required the characteristics i d (λ d , λ q , θ i ) and i q (λ d , λ q , θ i ) to be in the inverse form; the currents are functions of fluxes and positions.

Inverse Solution of Currents versus Flux Linkages
Equation (11) is used to represent the error between required fluxes (λ do and λ qo ) and actual fluxes, λ d (i d , i q , θ i ) and λ q (i d , i q , θ i ), that correspond to currents i d and i q . The currents i d and i q are obtained by minimizing the objective function F obj using an iterative  Figure 7 shows the flowchart of this process. First, the desired or commanded fluxes (λ do and λ qo ) are inputted. Then, the currents (i d and i q ) are changed with small steps (∆i d = 0.1 A and ∆i q = 0.1 A). for each value of currents (i d and i q ) within the iterative process, the fluxes are estimated from the flux maps of Figure 6. Finally, the objective function is estimated for obtained current and flux vectors, and then the currents (i d and i q ) that corresponds to minimum flux error (minimum F obj ) are output as the optimal solution of currents. Figure 8 shows the obtained optimal i d and i q currents after the Inverse solution of currents versus flux linkages.

Production of MATLAB Simulink Model
The proposed machine model that is given in Figure 2 consists of three main parts: first, the torque estimation Te (id, iq, θi), which is a straightforward task, as the obtained data from FEA can be stored in the form of a 3D LUT; second, the representation of iron estimated from the flux maps of Figure 6. Finally, the objective function is estimated for obtained current and flux vectors, and then the currents (id and iq) that corresponds to minimum flux error (minimum Fobj) are output as the optimal solution of currents. Figure 8 shows the obtained optimal id and iq currents after the Inverse solution of currents versus flux linkages.
start Calculate and save flux linkages (λd(id,iq,θ), λq(id,iq,θ)) from flux maps (Figure 6a, b) The obtained id and iq currents of Figure 8 are inputted to flux maps (Figure 6a,b). The resultant fluxes are compared to their corresponding flux values. The flux errors are shown in Figure 9. As noted, the maximum error is very small and can be ignored. The small error reflects the accuracy of the inverse solution and hence the accuracy of the proposed model. Calculate and save flux linkages (λd(id,iq,θ), λq(id,iq,θ)) from flux maps (Figure 6a, b) The obtained id and iq currents of Figure 8 are inputted to flux maps (Figure 6a,b). The resultant fluxes are compared to their corresponding flux values. The flux errors are shown in Figure 9. As noted, the maximum error is very small and can be ignored. The small error reflects the accuracy of the inverse solution and hence the accuracy of the proposed model. The obtained i d and i q currents of Figure 8 are inputted to flux maps (Figure 6a,b). The resultant fluxes are compared to their corresponding flux values. The flux errors are shown in Figure 9. As noted, the maximum error is very small and can be ignored. The small error reflects the accuracy of the inverse solution and hence the accuracy of the proposed model.

Estimation of Iron Losses
The iron losses per unit volume can be estimated by the well-known Bertotti formula, as a function of frequency (f) and maximum flux density (Bm), as given by Equation (12) [34,35]. kh, kc, and ke are the coefficients of hysteresis, eddy current, and

Estimation of Iron Losses
The iron losses per unit volume can be estimated by the well-known Bertotti formula, as a function of frequency (f ) and maximum flux density (B m ), as given by Equation (12) [34,35]. k h , k c , and k e are the coefficients of hysteresis, eddy current, and additional losses, respectively. k c is calculated as a function of material conductivity (σ) and lamination thickness (k d ), as given by Equation (13).
Equation (14) is a reformulation of Equation (12). The flux densities are converted to their corresponding flux-linkages based on machine dimensions and winding information. The equivalent iron loss resistance (R c ) is given by Equation (16).
k c = π 2 σk 2 d /6 (13) where V t and V y , h t and h y , and A t and A y are the volumes, heights, and areas for the stator tooth and yoke, respectively. The equivalent circuit of an IPMSM representing iron losses P fe and copper losses P cu is demonstrated in the rotor reference frame as shown in Figure 10. The iron and copper losses are represented by the resistances R c and R s , respectively. i da and i qa are the magnetizing currents. They can be calculated using Equation (17). The torque expression in terms of i da and i qa is given by Equation (18).

Model Verification
The proposed model, shown in Figure 2, was implemented in the MATLAB/Simulink environment. It is validated in two ways. The first was a comparison with the FEA model, and the second was by the experimental measurements on a physical 12/10 IPMSM prototype. Moreover, the experimental verifications were achieved in both the generator and motor actions.

Compared to FEA Model
To check the accuracy of proposed model, first, the Simulink model was run at a

Model Verification
The proposed model, shown in Figure 2, was implemented in the MATLAB/Simulink environment. It is validated in two ways. The first was a comparison with the FEA model, and the second was by the experimental measurements on a physical 12/10 IPMSM prototype. Moreover, the experimental verifications were achieved in both the generator and motor actions.

Compared to FEA Model
To check the accuracy of proposed model, first, the Simulink model was run at a specific operating point based on the field-oriented control (FOC) algorithm, and the resulting current (i d , i q ) waveforms were extracted. These currents (i d , i q ) were then injected into the FEA model to estimate the fluxes and electromagnetic torque.
The comparison was achieved at speed of 1500 r/min and stator current of 10 A. The waveforms of current (i d , i q ) are as shown in Figure 11a,b. The d-and q-axes flux linkages were estimated by the FEA model and then compared to corresponding values calculated within the MATLAB Simulink model. The comparison is given in Figure 11c,d. It was observed that the FEA-predicted d-and q-axes fluxed as well as the torque, as seen in Figure 11e, and could be completely reproduced by the MATLAB Simulink model. Note that the Simulink model required a very short time compared to FEA to produce the same results.
The proposed model, shown in Figure 2, was implemented in the MATLAB/Simulink environment. It is validated in two ways. The first was a comparison with the FEA model, and the second was by the experimental measurements on a physical 12/10 IPMSM prototype. Moreover, the experimental verifications were achieved in both the generator and motor actions.

Compared to FEA Model
To check the accuracy of proposed model, first, the Simulink model was run at a specific operating point based on the field-oriented control (FOC) algorithm, and the resulting current (id, iq) waveforms were extracted. These currents (id, iq) were then injected into the FEA model to estimate the fluxes and electromagnetic torque.
The comparison was achieved at speed of 1500 r/min and stator current of 10 A. The waveforms of current (id, iq) are as shown in Figure 11a,b. The d-and q-axes flux linkages were estimated by the FEA model and then compared to corresponding values calculated within the MATLAB Simulink model. The comparison is given in Figure  11c,d. It was observed that the FEA-predicted d-and q-axes fluxed as well as the torque, as seen in Figure 11e, and could be completely reproduced by the MATLAB Simulink model. Note that the Simulink model required a very short time compared to FEA to produce the same results.

Experimental Verification as a Generator
To verify the effectiveness of the proposed Simulink model, an experimental platform was set up, as shown in Figure 12. The drive system included a three phase inverter, an IGBT gate-driver, current transducers (LEM LA100-P), voltage transducers (LEM LV 25-P), a torque sensor (Lorenz DR 2112), an incremental encoder (1024 PPR), a data acquisition board (NI 6229 DAQ), an electromagnetic brake, and a digital signal processor (C6713 DSP). The DSP was connected to an FPGA for better timing. The data were collected online with the help of a daughter card using MATLAB.
The IPMSM was operated as a generator. It was driven by a dc motor. First, Figure  13 shows the experimental results under no load at 1708 r/min. The experimental results showed three phases' currents, the three phases' voltages, and electromagnetic torque.

Experimental Verification as a Generator
To verify the effectiveness of the proposed Simulink model, an experimental platform was set up, as shown in Figure 12. The drive system included a three phase inverter, an IGBT gate-driver, current transducers (LEM LA100-P), voltage transducers (LEM LV 25-P), a torque sensor (Lorenz DR 2112), an incremental encoder (1024 PPR), a data acquisition board (NI 6229 DAQ), an electromagnetic brake, and a digital signal processor (C6713 DSP).
The DSP was connected to an FPGA for better timing. The data were collected online with the help of a daughter card using MATLAB.
processor (C6713 DSP). The DSP was connected to an FPGA for better timing. The data were collected online with the help of a daughter card using MATLAB.
The IPMSM was operated as a generator. It was driven by a dc motor. First, Figure  13 shows the experimental results under no load at 1708 r/min. The experimental results showed three phases' currents, the three phases' voltages, and electromagnetic torque. The average value of no-load torque, as seen in Figure 13c, was 0.39 Nm, which basically represents iron and friction losses. The torque ripple was very noticeable. Figure 13d shows a comparison between the measured phase voltage and the phase voltage obtained from FEA. As observed, a perfect match is seen in Figure 13d that reflects the accuracy of the FEA model and hence the accuracy of the proposed model.  The IPMSM was operated as a generator. It was driven by a dc motor. First, Figure 13 shows the experimental results under no load at 1708 r/min. The experimental results showed three phases' currents, the three phases' voltages, and electromagnetic torque. The average value of no-load torque, as seen in Figure 13c, was 0.39 Nm, which basically represents iron and friction losses. The torque ripple was very noticeable. Figure 13d shows a comparison between the measured phase voltage and the phase voltage obtained from FEA. As observed, a perfect match is seen in Figure 13d that reflects the accuracy of the FEA model and hence the accuracy of the proposed model. connected to a resistive load. Figure 14d shows a comparison between the measured phase voltage and the phase voltage obtained from FEA. As observed, a perfect match was also seen.  Second, Figure 14a-c show the experimental results under loading conditions as a generator. The IMPSM was driven at a speed of 854 r/min, and the terminals were connected to a resistive load. Figure 14d shows a comparison between the measured phase voltage and the phase voltage obtained from FEA. As observed, a perfect match was also seen.  Figure 15a shows a comparison between the measured torque (by a torque sensor) and the predicted FEA torque under loading conditions of Figure 14. As seen, the torque sensor was not fast enough to capture all the torque details. The full torque details, including torque ripple, was well displayed by the FEA model. Moreover, Figure 15b gives the average torque values. The difference between the measured and FEA came  Figure 15a shows a comparison between the measured torque (by a torque sensor) and the predicted FEA torque under loading conditions of Figure 14. As seen, the torque sensor was not fast enough to capture all the torque details. The full torque details, including torque ripple, was well displayed by the FEA model. Moreover, Figure 15b gives the average torque values. The difference between the measured and FEA came from the fact that FEA does not consider iron loss effect. The FEA provides a magnetostatic solution. The iron loss effect was not considered in the FEA model as it is usually evaluated in post processes of most commercial FEA tools. On the other hand, the iron loss was included within the proposed model. Hence, very good agreement was observed between the measured torque signal and the simulated torque signal, as seen in Figure 15b. from the fact that FEA does not consider iron loss effect. The FEA provides a magnetostatic solution. The iron loss effect was not considered in the FEA model as it is usually evaluated in post processes of most commercial FEA tools. On the other hand, the iron loss was included within the proposed model. Hence, very good agreement was observed between the measured torque signal and the simulated torque signal, as seen in Figure 15b. Table 3 shows the analysis details for torque data. As illustrated, the torque error for the FEA model was 18.8%. Hence, the FEA model was not accurate enough to represent the real machine, as it neglects iron loss effect. On the contrary, the proposed Simulink model provided a torque error of 2.13%, which may have been a result of stray losses, variation of dc voltage, measurement errors (signal errors and quantization errors), tolerances in machine dimensions, and degradation of material properties. Hence, the proposed model succeeded in representing the real machine accurately.     Table 3 shows the analysis details for torque data. As illustrated, the torque error for the FEA model was 18.8%. Hence, the FEA model was not accurate enough to represent the real machine, as it neglects iron loss effect. On the contrary, the proposed Simulink model provided a torque error of 2.13%, which may have been a result of stray losses, variation of dc voltage, measurement errors (signal errors and quantization errors), tolerances in machine dimensions, and degradation of material properties. Hence, the proposed model succeeded in representing the real machine accurately. Table 3. Summary of torque results.

Experimental Verification as a Motor
Finally, experimental tests for the IPMSM in motoring action were achieved based on the FOC algorithm. The dc link voltage was set to 200 V. First, the FOC algorithm was tested to ensure proper tracking operation. Figure 1 shows the performance of the FOC algorithm under sudden changes of reference currents (i d-ref and i q-ref ). As seen in Figure 16a,b, i q-ref changed as a triangular wave at a frequency of 20 Hz. The motor speed was 486 r/min. As seen in Figure 16c,d, both i d-ref and i q-ref changed as square waves at a frequency of 10 Hz. The motor speed was 488 r/min. As noted, for both cases, the FOC algorithm worked efficiently, as the motor could track the reference d-and q-axes currents properly. Moreover, the three phases currents were symmetric and sinusoidal.   Table 4 gives the torque analysis data. It compares the measured and the simulated torques. The torque error was 3.26% which may be a result of measurement errors. Moreover, the mechanical coupling between the IPMSM and loading magnetic brake could be a primary reason for torque error, as improper mechanical coupling could affect friction losses. As a result, the proposed model is accurate enough to represent real machine performance. Recently, many researchers have become interested in MTPA and maximum efficiency per ampere (MEPA) control for IPMSM, as they are very important, especially for high power machines such as in EV applications. The basic idea for MTPA and MEPA is to define the current angle that corresponds to the desired criteria (maximum torque or maximum efficiency). The proposed model can easily be used for such purposes. The obtained current angles can be further employed in real implementations to achieved MTPA and MEPA. Figure 18 shows a comparison of the measured and simulated torque versus current angle. The current angle was changed from 0 • to 50 • in steps of 5 • . The maximum current was set to 3 A. The motor was running in an open loop current control. The motor speed was adjusted by the loading dc brake. The test was done at 500 r/min. As noted, due to the reluctance torque of IPMSM, the maximum torque did not occur at a zero-current angle, such as in surface-mounted PMSM. As seen, the maximum torque occurred around a 5 • current angle. Moreover, very good agreement was observed between measurement and simulation, as seen in Figure 18a. The maximum torque error over 11 points of comparison was 0.12 Nm at a current angle of 0 • , as illustrated by Figure 18b. It corresponded to an error of 4.1%, as seen in Figure 18c. On the other hand, the average torque error over 11 points of comparison was 0.0449 Nm, which corresponded to a percentage of 1.7422%, which is a small error and could be a result of measurement errors, assumptions in the FEA model, variations of dc voltage, or effects of temperature on machine parameters such as magnets and resistance. current was set to 3 A. The motor was running in an open loop current control. The motor speed was adjusted by the loading dc brake. The test was done at 500 r/min. As noted, due to the reluctance torque of IPMSM, the maximum torque did not occur at a zero-current angle, such as in surface-mounted PMSM. As seen, the maximum torque occurred around a 5˚ current angle. Moreover, very good agreement was observed between measurement and simulation, as seen in Figure 18a. The maximum torque error over 11 points of comparison was 0.12 Nm at a current angle of 0°, as illustrated by Figure 18b. It corresponded to an error of 4.1%, as seen in Figure 18c. On the other hand, the average torque error over 11 points of comparison was 0.0449 Nm, which corresponded to a percentage of 1.7422%, which is a small error and could be a result of measurement errors, assumptions in the FEA model, variations of dc voltage, or effects of temperature on machine parameters such as magnets and resistance.
The model accuracy will remain similar for wider speed ranges as long as the control algorithm (FOC) provides good tracking behavior. If the FOC fails tracking, this will not affect the model accuracy, which means that the control algorithm (FOC) itself needs improvements, not the proposed machine model.  Figure 19 shows a flowchart that explains step by step the proposed modeling method.

Summary of Proposed Method
Enter full dimensions of machine geometry, Enter material properties for magnets, steel, and windings, Define electric circuits and their configurations.
Construct the finite element model for tested machine.
Check the accuracy of finite element model by comparison with the measured voltages as in Figures 13d and 14d.

Accuracy is accepted?
yes Calibrate finite element model by the measurement of PM flux linkage (Section 3.1.1)

No
Generate the full magnetic characteristics for tested machine using FEA such as in Figure 6.
Inverse solution of currents versus flux-linkages (section 3.2.1) Build the MATLAB Simulink model of Figure 2 involving iron losses of Section 3.2.2 Validate model accuracy experimentally in motoring action such as in Figures 18 and 19 Accuracy is accepted? The model accuracy will remain similar for wider speed ranges as long as the control algorithm (FOC) provides good tracking behavior. If the FOC fails tracking, this will not affect the model accuracy, which means that the control algorithm (FOC) itself needs improvements, not the proposed machine model. Figure 19 shows a flowchart that explains step by step the proposed modeling method. First, the machine dimensions, material properties, and configurations of current circuits are defined. Second, the FE model of the machine is constructed. The FE model must be precisely configured to represent real machine behavior. Any mismatches in input data to the FE model will lead to undesirable and huge errors. Hence, the experimental verification of the FE model is of great importance. Third, the FE model is validated by comparing the no-load and loading voltages of the machine as a generator. If the accuracy is accepted and hence the PM flux linkage is inputted precisely, then one can proceed to calculate the full magnetic characteristics of machine. on the contrary, if the accuracy is not satisfactory, then measurement of the PM flux linkage must be done precisely after an FFT of phase voltage, as described in Section 3.1.1. After that, the production of the full magnetic characteristics of the machine can be achieved. Fourth, the development of the MATLAB Simulink model is done using the FEA data after precise calibrations. The MATLAB Simulink model requires an inverse solution of currents versus flux linkages, as discussed in Section 3.2.1. It also involves the estimation of iron losses and its involvement within the MATLAB Simulink model, as described in Section 3.2.2. The final step is to validate the MATLAB Simulink model experimental when the machine is operated as a motor.

Summary of Proposed Method
(c) Figure 18. Comparison between the measured and simulated torque vs. current angle: (a) measured and simulated torques vs. current angle; (b) torque error in Nm; (c) percentage torque error. Figure 19 shows a flowchart that explains step by step the proposed modeling method.

Summary of Proposed Method
Enter full dimensions of machine geometry, Enter material properties for magnets, steel, and windings, Define electric circuits and their configurations.
Construct the finite element model for tested machine.
Check the accuracy of finite element model by comparison with the measured voltages as in Figures 13d and 14d.

Accuracy is accepted?
yes Calibrate finite element model by the measurement of PM flux linkage (Section 3.1.1)

No
Generate the full magnetic characteristics for tested machine using FEA such as in Figure 6.
Inverse solution of currents versus flux-linkages (section 3.2.1) Build the MATLAB Simulink model of Figure 2 involving iron losses of Section 3.2.2 Validate model accuracy experimentally in motoring action such as in Figures 18 and 19 Accuracy is accepted?  Figure 19. Flowchart of the proposed modeling method.
First, the machine dimensions, material properties, and configurations of current circuits are defined. Second, the FE model of the machine is constructed. The FE model must be precisely configured to represent real machine behavior. Any mismatches in input data to the FE model will lead to undesirable and huge errors. Hence, the experimental verification of the FE model is of great importance. Third, the FE model is validated by comparing the no-load and loading voltages of the machine as a generator. If the accuracy is accepted and hence the PM flux linkage is inputted precisely, then one can proceed to calculate the full magnetic characteristics of machine. on the contrary, if the accuracy is not satisfactory, then measurement of the PM flux linkage must be done precisely after an FFT of phase voltage, as described in Section 3.1.1. After that, the production of the full magnetic characteristics of the machine can be achieved. Fourth, the development of the MATLAB Simulink model is done using the FEA data after precise calibrations. The MATLAB Simulink model requires an inverse solution of currents versus flux linkages, as discussed in Section 3.2.1. It also involves the estimation Figure 19. Flowchart of the proposed modeling method.

Conclusions
This paper develops an accurate modeling method for PMSMs. First, the proposed model exhibits high fidelity in predicting real machine behavior while providing a computationally efficient solution. It presents a powerful solution for detailed analysis of all kinds of PMSMs for further performance analysis or for developing accurate control algorithms. It also helps to analyze machine losses and estimate real machine efficiency. Accurate loss modeling and efficiency estimation are key interests for current and future research related to several industrial applications such as electric vehicles. Second, the proposed model considers accurately iron losses effects as well as the effects of magnetic saturation, spatial harmonics, and mutual coupling. Third, it not only guaranties the production of a high accuracy model but also involves the manufacturing effects by the calibration of the FEA model based on experimental measurements. As seen by the results, the good match between simulation and experimental measurements reflects the high model accuracy. The average percentage of torque error between measurement and simulation is 1.7422%, which is a small and acceptable error.

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

Nomenclature Symbol
Definition Unit A t , A y The areas of stator pole/tooth and yoke m 2 B The combined rotor and load viscous friction coefficient kg·m 2 B m The peak value of magnetic flux density Wb/m 2 f The frequency of back EMF voltage Hz F obj The objective function h t , h y The heights of stator pole/tooth and yoke m i d , i q The d-and q-axes components of stator current A i da , i qa The d-and q-axes components of magnetizing currents A i d-fe , i q-fe The d-and q-axes components of currents that represents iron losses A i d-ref , i q-ref The d-and q-axes components of reference current A J The combined rotor and load inertia coefficient kg·m 2 ·s k c The coefficients of eddy current losses W/m 3 k d The lamination thickness m k e The coefficients of additional losses W/m 3 k h The coefficients of hysteresis losses W/m 3 L d , L q The d-and q-axes components of inductance H p The number of pole pairs P cu The copper losses W P fe The iron losses W R c The iron loss resistance Ω R s The phase resistance Ω T e The total electromagnetic torque N·m T L The load torque N·m v d , v q The d-and q-axes components of phase voltage V V t , V y The volumes of stator pole/tooth and yoke m 3 V ph_peak The peak value of phase voltage V