Novel ABC—Variable Based Closed-Form Analytical Model of Surface Mounted Permanent Magnet AC Machines With Sinusoidal and Non-Sinusoidal Back-EMF

This work presents a novel abc-based model applicable to surface-mounted permanent magnet AC (SM-PMAC) machines with sinusoidal and non-sinusoidal back-electromotive force (back-emf). It is capable of predicting the electromagnetic performance metrics such as torque waveforms, machine inductances, flux linkages and back-emf. The closed form expressions of the model, which can be evaluated with a high computational efficiency, are derived from basic geometric and winding parameters. Validation of the model is carried out numerically and experimentally with a very good match in results. Finally, the computational efficiency of the model is highlighted by considering a multi-objective evolutionary optimization design of SM-PMAC machine with a relatively large number of design parameters, where results are presented and discussed.

models. The accuracy of predictions made from FEA was 23 shown to be comparable to experimental measurements [1], 24 [2], [3], [4]. However, this method suffers from computational 25 inefficiency and it might be impractical to use in the initial 26 The associate editor coordinating the review of this manuscript and approving it for publication was Yilun Shang. stages of the design when many iterations are required to 27 study the variation of different parameters on the machine 28 performance [5], [6]. 29 Analytical models on the other hand can be convenient 30 for this purpose as they are much faster to evaluate and 31 their accuracy has been steadily improving [ [12]. For instance, in [7], a high-speed permanent 33 magnet synchronous machine (PMSM) was designed using 34 an analytical model that in addition to the electromagnetic 35 performance; it could predict important thermal and mechan-36 ical metrics. A machine with a very good match between 37 predicted and experimental measurements was designed and 38 tested. In [8], an analytical model of cogging torque in 39 PMSMs was developed. This model was capable of cap-40 turing the airgap flux density with results comparable to 41 FEA considering eccentricity in the rotor or defects in the 42 magnets. It took just 15 seconds to solve the analytical model 43 in contrast to one hour in FEA. 44 Sub-domain models were also applied for analyzing elec-45 tric machines in a number of research work [13], [14], [15], 46 [16], [17]. In these methods, the machine is divided into subregions where Maxwell's equations are developed in each • The analytical model presented is valid for PMAC 92 machines with sinusoidal and non-sinusoidal back-emf. analysis, and therefore can be easily applied to analyze 101 a given PMAC machine. 102 The organization of this paper is as follows. In Section II 103 the model is presented and derived in terms of basic machine 104 geometrical and winding parameters. In Section III the model 105 is validated using numerical and experimental measurements 106 on a lab-scale prototype. In Section IV, the model is used in 107 an evolutionary multi-objective optimization environment to 108 demonstrate its effectiveness from a computational perspec-109 tive. In Section V, the conclusion and future work drawn from 110 this work are presented.

112
The proposed model is derived and explained in this section. 113 In this derivation, the following assumptions were made. 114 First, the magneto-motive force (MMF) drop across the steel 115 is negligible. Second, only the radial component of the airgap 116 flux density is considered, the tangential and axial compo-117 nents are neglected due to their negligible magnitude. Third, 118 the three-phase stator currents are assumed to be balanced and 119 harmonic free. Lastly, tooth-tips in the stator teeth are ignored 120 and the teeth are assumed to have a rectangular geometry.

122
A developed diagram of a surface-mounted PMAC machine, 123 which shows a portion of it redrawn in a linear fashion, 124 is shown in Fig. 1. The spatial mechanical position of the 125 machine φ sm is defined with respect to the stator reference 126 axis located at the center of Stator Tooth 1 (ST1). Based on 127 this, the rotor position θ rm is defined as the displacement with 128 respect to the stator reference axis. Note that blue sections 129 indicate magnets whose flux points radially outward from 130 rotor to stator, while red sections are the magnets with flux 131 pointing radially inward from stator to rotor. Each magnet 132 spans a mechanical angle θ pm given by where P r is the number of magnet poles and α pm is a fraction 135 between 0 and 1, controlling the circumferential span of the 136 magnets. Some additional crucial parameters to the design 137 and optimization process that are labeled in Fig. 1, are the 138 depths of the stator back-iron d sb and teeth d tb , the depth of 139 the rotor back-iron d rb , the depth of the magnets d m , airgap 140 length g, the depth of the inert region d i between shaft and 141 rotor and the radius of rotor shaft r rs . Each stator tooth has a 142 uniform width w tb . Using these variables, a number of param-143 eters such as mass and volume of different components can be 144 derived [27], [28]. A more detailed geometrical description of 145 the machine can be found in [18].

147
The stator winding distribution for each phase can be 148 expressed using the winding function, which is a continuous 149 function giving the number of turns around each stator tooth. 150 VOLUME 10, 2022 The a-phase winding function can be expressed as where w as is a-phase continuous winding function, k is an working harmonic is such that an abc sequence is established, 163 parameter γ is added, which is set to 1 or −1 according to the 164 condition explained in [27]. Note that the series summation 165 over k is not shown in (2) for clarity.

167
The a-phase current can be expressed as where I s is the rms current and φ i is the current vector position tively. A closed-form expression for (4) is presented in [27] 180 and [28] for concentrated winding and sinusoidally dis-181 tributed winding distributions; respectively.

182
The MMF due to the rotor magnets can be expressed using 183 an odd function as where n is the magnet MMF odd positive harmonic number, 186 and f pm,n is the harmonic component given by with µ 0 is the vacuum permeability, µ r is the magnet relative 189 permeability and B r is the residual flux density. Note that the 190 summation sign to represent the series is omitted for clarity 191 here and in the following equations. Therefore, as presented in [29] and [30], the maximum 202 possible airgap depth to which the flux can reach in a stator 203 slot (SS), d fmx , depends on the width of the slot w so , the depth 204 of the magnets d m and their relative permeability µ r , and the 205 length of the airgap g. It can be expressed as Hence, and as demonstrated in Fig. 2, it is possible to define 208 an effective tooth depth d * tb representing the maximum depth 209 of the slot the flux travels before emerging to the teeth given 210 by Based on the previous analysis, the airgap permeance vari-213 ation shown in Fig. 2 can be expressed using the Permeance 214 function as where j is an even positive integer, S s is the number of stator 217 teeth, and p dc and p ac j are the dc and ac harmonic components, 218 respectively; given by the following two expressions 219 p dc = P g min + P g max − P g min α t (10) 220 where α t is a fraction, theoretically in the range of 0 and 1, 222 controlling a stator tooth angular span θ t . This angle along 223 with P g,min and P g,max are defined as where r st is the radius to a stator tooth, and r rb is the radius 229 to the outer region of the rotor back-iron.

230
Based on (4), (5) and (9), and using a radial field analy-231 sis [27], [28], the airgap flux density spatial and temporal 232 distribution at any radius r between the outer rotor back-iron 233 and stator teeth tips can be predicted using Considering a balanced three-phase system, the abc induc-238 tances can be obtained utilizing the following expression where l is the stack length of the machine, w xs is the x-phase were presented previously but assumed that the airgap is 246 spatially uniform [19]. In this work, to increase the accuracy; 247 the spatial variation is also considered. Using Fourier series, 248 this function can be expressed as an even function given by where h is an even positive harmonic number, and g dc and g h 251 are the dc and ac harmonics components, respectively; given 252 by [31] The inverse of (17) is defined as e v , which can be shown to 256 be equal to with e 1 and e 2 are given by where L ls is the self-leakage inductance considering both slot-274 windings and end-windings leakages as detailed in [18], L dc 275 and L ac are the dc and ac-components of the self-inductance 276 given by the following two expressions where N pc is the number of turns per coil and w s is obtained 280 by squaring the coefficient w k for each harmonic and then 281 summing them as follows Similar expressions can be obtained for b-phase and 284 c-phase with the appropriate phase shifts. In the same man-285 ner, applying (16) for the mutual inductances between any 286 two phases, expressions for the mutual inductances can be 287 VOLUME 10, 2022 obtained. For instance, the mutual inductance between the a-288 phase and b-phase is given by where L mdc is the dc component of the mutual-inductances   Substituting the expressions for currents (3), induc-348 tances (16), and magnet flux linkages (32) into (40), and 349 from (41); the torque expression can be obtained. It can be 350 shown that the torque expression can be segregated into two 351 terms: a constant independent of rotor position term T L which 352 stems from the first term on the right-hand side of (40), and 353 coefficient term that is multiplied by terms that depend on the 354 rotor position resulting from the second term in (40): where T L , T e0 , T e1 , and T e2 are given by 357 T L = 9 4 I 2 s P r L ac sin (2φ i ) To get more insight on (42), Fig. 3 shows the predicted

389
Details on the machine geometry, magnet, winding and 390 applied current are given in Table 1, where I rated is the rated 391 line current, V rated is the rated line-line voltage, P rated is the 392 rated output power, T rated is the rated torque, and ω rated is the 393 rated speed. The winding function harmonic coefficient for 394 this machine is given in the first row of Table 2.    shows that such assumption can introduce some error in the 427 predictions.  analytical models are found to be 0.91% and 0.75% for op1, 435 0.90% and 0.74% for op2, and 0.93% and 0.76% for op3, 436 respectively. Nevertheless, a good match in results can be 437 seen.

438
It is noted that the predicted error increases as the excita-439 tion level increases, and that the predicted analytical torque 440 is higher than that of the FEA This error stems from the 441 assumption made in the analytical model which ignores the 442 tangential component in the airgap flux density and assumes 443 all of the flux density in the airgap is radial. As the excita-444 tion level keeps increasing, the tangential component in the 445 airgap increases. The analytical model assumes that all flux 446 density in the airgap is radial and uses that to calculate the 447 torque, which results in overestimation when compared to 448 FEA torque.

449
The airgap flux density under full rated torque is shown in 450 Fig. 6. It is observed that slotting effects are captured with 451 good accuracy.

453
Experimental verification was carried out using the test-setup 454 shown in Fig. 7 and with the parameters listed in Table 1. This 455 setup consisted of a Y-configured PMSM coupled to a DC 456 machine. The DC machine acts as generator and the output of 457 which is connected to a resistive load. The PMSM is powered 458 from an electronic variable frequency drive which generates 459 the three-phase PWM signal in an open-loop mode. A torque 460 sensor, rated at 100 Nm, 10000 rpm and with sensitivity of 461 50 mV/Nm, mounted on the shaft is used for torque measure-462 ment. The machine stator inductance is measured using an 463 LCR meter. Table 4 shows a comparison of the average torque val-465 ues obtained analytically and experimentally for the three 466 aforementioned operating points. Note that due to equip-467 ment limitations, torque ripple could not be measured. 468 Good agreement in results can be seen. The experimentally 469 measured a-phase self-inductance was equal to 2.20 mH 470 which is close to the analytical model predicted value of 471 2.26 mH shown in Table 3. The no-load back-emf wave-472 forms at the rated conditions obtained from the analyti-473 cal model, FEA and experimental measurement are shown 474 in Fig. 8 where it is noted that the back-emf is sinu-475 soidal for these conditions. Good agreement in results 476 is seen.   function and therefore the amplitude of the winding function 483 harmonic coefficients, and the slot and poles combination. 484 Based on this, an FEA model with magnet poles equal to 485 30 and with α t and α pm set 0.5 and 0.7, respectively; was 486 built. The rest of the machine parameters are unchanged as 487 listed in Table 1 while the winding configuration applied has 488 the winding function harmonics coefficient listed in row 2 of 489  Fig. 9 shows the back-emf obtained from the developed 493 model and FEA. As shown, the back-emf is non-sinusoidal 494 due to the existence of significant higher order harmonics 495 in (33). The analytically predicted average torque, torque 496 ripple percentage and a-phase self-inductance at θ rm equals 497 to zero were calculated at 15.55 Nm, 0.8%, and 8.78 mH, 498 respectively. The same values obtained from FEA were found 499 to be 14.75 Nm, 1.1%, and 8.12 mH; respectively. The 500 increased error in the average torque prediction compared to 501 the results in sub-sections A and B is due to the relatively 502 more magnet-magnet leakage present in Vernier machine 503 compared to FSCW-PMSM.

505
The developed analytical model was integrated into an evo-506 lutionary optimization design algorithm with the purpose 507 of designing a 1.86 kW PMSM with a 5:1 constant power 508 speed range (CPSR) [27]. In this design, three torque-speed 509 operating points were modeled to capture the flux weakening 510 VOLUME 10, 2022 region [27]. A total of 20 design parameters/degrees of free-511 dom were used [27]. Ten of these parameters were geomet-512 rical and winding parameters, four were used to select the 513 material of the stator, rotor, magnets, and conductors; and the 514 last six were the stator current rms value and current angle 515 for each operating point [27]. A total of 18 hard constraints 516 were applied and two optimization objectives were set to 517 be minimized which were the total machine electromagnetic 518 mass and total weighted loss computed from the three points.

519
The total weighted loss includes the winding conduction loss, 520 core loss, magnet loss and conduction loss in the inverter. topology is adopted as described in [27]. More information 534 about this design procedure can be found in [27].

535
The optimization was run for a GA population of 500 and 536 generation of 500. The resulting trade-off curve between the 537 defined objectives, better known as the Pareto-optimal front, 538 is shown in Fig. 10, where the x-axis represents the total 539 electromagnetic mass of the machine and y-axis gives the 540 total weighted power loss. Each point in this curve is a sep-541 arate machine design that satisfied all imposed constraints. conduction and core loss happening in the stator region, but 550 at the expense of heavier machines.

551
A total of 147 machines were designed as shown in Fig. 10, 552 where the optimization took around one hour on an i-5 desk-553 top PC with 8 GB-ram, 3.0 GHz processor, and 8 cores, which 554 involved 1.2 million evaluations of the analytical model 555 script. Similar design procedures using numerical methods 556 with much smaller space require much more computational 557 time [5], [6], [37].

558
To assert the computational effectiveness of the developed 559 model, it was run using 300 winding function harmonics k 560 in (2), 300 magnet MMF harmonics n in (5), and 300 airgap 561 permeance harmonics j in (9). Moreover, 1200 points were 562 considered, in the span of zero to 2π, for φ sm . The analytical 563 script was run for one rotor position corresponding to θ rm 564 equal to zero. To make the comparison as fair as possible, 565 two magneto-static linear FEA studies were conducted for the 566 same rotor position: the first incorporated adaptive meshing 567 with three pass/adaptions; while the second did not incorpo-568 rate adaptive meshing and was relying on the initial mesh. 569 By running the studies on an i-5 desktop PC with 3.0 GHz and 570 8 processing cores, it was found that the analytical script took 571 0.025 seconds to obtain the electromagnetic field solutions, 572 compared to 15 seconds and 9 seconds for the first and 573 second FEA studies, respectively. It should be noted that this 574 time does not include the time needed to build the machine 575 geometry or its winding in the analytical or FEA software.

577
This work presents a novel abc-variable based analytical 578 model for surface-mounted PMAC machines which can be 579 utilized to predict the electromagnetic machine performance 580 with a sinusoidal or non-sinusoidal back-emf. Closed-form 581 expressions for important parameters such as inductances, 582 flux linkages, and torque were presented. Validation of the 583 model was done using FEA and experimental tests on a lab-584 prototype, where a good match in results was demonstrated 585 for sinusoidal and non-sinusoidal back-emf. In particular, 586 emf, the maximum error in the average torque between the  Finally, the inclusion of cogging torque prediction can be 601 considered.