Exact Two-Dimensional Analytical Calculations for Magnetic Field , Electromagnetic Torque , UMF , Back-EMF , and Inductance of Outer Rotor Surface Inset Permanent Magnet Machines

This paper presents a two-dimensional analytical model of outer rotor permanent magnet machines equipped with surface inset permanent magnets. To obtain the analytical model, the whole model is divided into the sub-domains, according to the magnetic properties and geometries. Maxwell equations in each sub-domain are expressed and analytically solved. By using the boundary/interface conditions between adjacent sub-regions, integral coefficients in the general solutions are obtained. At the end, the analytically calculated results of the air-gap magnetic flux density, electromagnetic torque, unbalanced magnetic force (UMF), back-electromotive force (EMF) and inductances are verified by comparing them with those obtained from finite element method (FEM). One of the merits of this method in comparison with the numerical model is the capability of rapid calculation with the highest precision, which made it suitable for optimization problems.


Introduction
The existence of different types of PM brushless machines (PMBLM) made them applicable for a wide range of applications.PMBLMs have superiorities in comparison with other rivals like induction machines or reluctance synchronous machines due to higher efficiency, high torque per volume, lower torque ripple, lower vibration, and lower acoustic noise.
PMBLM can be categorized in terms of various criteria such as the topology of PMs, the relative position of the rotor and stator, the slotted or slotless stator structure, etc.
Various PM topologies such as surface-mounted, surface-inset, and interior are used where each has its own advantages and disadvantages.Among these topologies, surface-inset can provide a compromise between the other two topologies.
Electric machines with single rotor and single stator can be either inner rotor or outer rotor.The outer rotor motors can develop more output torque than the inner ones for the same volume of the machine.Usually, inner rotor machines are used for applications, which need rapid acceleration and deceleration.Outer rotor machines usually are used for applications which need constant speed.Also, the mechanical robustness of the PMs in the outer rotor configuration is higher than the inner one.
In this paper because of the aforementioned advantages of the outer rotor machines and surface inset PMs, an exact two-dimensional electromagnetic model for this type of machines is extracted.
In the design procedure, the static model is normally considered.Numerous static models have been presented for electric machines, in which some of them are based on the analytical approaches , and the others are based on numerical methods [53][54][55].
The presented analytical model for electric machines are based on permeance model [1] or magnetic equivalent circuit, also known as (a.k.a.) 0-D analytic model [2], or resolution of the Maxwell's equations in 2-D plane (a.k.a.2-D analytic model)  or 3-D plane (a.k.a.3-D analytic model) [49,50].Also, other methods based on mapping techniques such as Schwarz-Christoffel have been used to extract the model of some machine with the analytical approach [51,52].The most accurate presented models among the analytical methods are 3-D analytic and 2-D analytic.The 2-D analytic method can be used instead of 3-D when the model is symmetric and has no skewing.Two-dimensional analytical models are not only capable of considering a high number of harmonics which elevate the precision of the models, but also have less computational time in comparison with the numerical methods and made them appropriate for optimal design problems.
Most of the developed 2-D or 3-D analytic model are assumed with the infinite permeability of the iron parts.New techniques to account for finite soft-magnetic permeability have been recently developed, i.e., the multi-layer model using the Cauchy's product theorem is presented in [38], and the subdomain technique by applying the superposition principle in both directions is proposed in [39][40][41][42][43][44][45][46] which can be used to calculate the core losses and saturation phenomena.
An overview of the analytical models in the Maxwell-Fourier method with a global or local saturation effect has been realized in [40].According to [48], Dubas' superposition technique [39,40] is very interesting since it enables the magnetic field calculation in the material of slotted geometries.This superposition technique has been implemented in radial-flux electrical machines with(out) PMs supplied by a direct or alternate current [44,45].
The presented technique in [39][40][41][42][43][44][45][46] is not only used to predict the magnetic field in all parts of the electrical machines, but also it is used to obtain a 2-D analytical model of the steady state heat transfer of the electrical machines by solving the heat equations [47].
The aim of this paper is to extract a 2-D analytical model of PM brushless outer rotor machines equipped with surface inset PMs.The model is used to analytically compute the electromagnetic torque, torque ripple, back-electromotive force (EMF), inductances and unbalanced magnetic forces (UMF).
Therefore, this paper is organized as follows.In Section 2 the procedure of extracting the 2-D model is explained.Section 3 is dedicated to the calculation of the electromagnetic quantities.In Section 4 the analytical results of the case study are presented and compared with those of the numerical method.In the final part this paper is concluded.

Assumptions
Figure 1 shows the topology of an outer rotor surface inset brushless permanent magnet machine.e) The edges of the slots and slot-openings have radial direction.f) The eddy current effect is neglected.

Dividing Region into Sub-regions
According to the shape and material characteristics, the whole domain is divided into a number of sub-domains.All the sub-domains are illustrated in Figure 1 and listed in Table 1 for a PMBLM with Q slots and p pole-pairs.
When the winding is single-layer or double-layer non-overlapping, as shown respectively in Figure 2a,b, each slot is a single subdomain; however, if the winding is two-layer overlapping, as shown in Figure 2c, each slot is divided into two subdomains.In order to make the problem solvable, below assumptions are made:

Dividing Region into Sub-Regions
According to the shape and material characteristics, the whole domain is divided into a number of sub-domains.All the sub-domains are illustrated in Figure 1 and listed in Table 1 for a PMBLM with Q slots and p pole-pairs.
When the winding is single-layer or double-layer non-overlapping, as shown respectively in Figure 2a,b, each slot is a single subdomain; however, if the winding is two-layer overlapping, as shown in Figure 2c, each slot is divided into two subdomains.

Extracting the Magnetic Model
In this part, for each sub-domain a partial differential equation is extracted based on Maxwell's equations.
Maxwell's equations in quasi-static form are as follows: where B is the magnetic flux density vector.Ampere's law represents the relation between the magnetic field intensity vector ( H ), external current density vector ( J ), and current density vector in media ( Ĵ ).In this investigation the current density vector in media is assumed to be negligible, i.e., ˆ0 = J .The relation between the magnetic flux density vector and the magnetic field intensity vector in permanent magnet media with linear demagnetizing curve is as follows: where M is the magnetization vector.
Substituting (3) in (2) yields The magnetic flux density vector can be represented as the curl of the magnetic vector potential ( A ): Using ( 4) and ( 5) the following expression is obtained: For each sub-domain, Equation (6) results in Poisson equations for the magnet and slot regions, and Laplace equations for the air-gap and slot-opening sub-domains, as represented below.
In 2-D polar coordinates, the magnetic vector potential and the current density vector have just a component along z , i.e., ( ) ( ) . Also, the magnetic flux density vector and the magnetization vector have radial and tangential components as below:

Extracting the Magnetic Model
In this part, for each sub-domain a partial differential equation is extracted based on Maxwell's equations.
Maxwell's equations in quasi-static form are as follows: where B is the magnetic flux density vector.Ampere's law represents the relation between the magnetic field intensity vector (H), external current density vector (J), and current density vector in media ( Ĵ).
In this investigation the current density vector in media is assumed to be negligible, i.e., Ĵ = 0.The relation between the magnetic flux density vector and the magnetic field intensity vector in permanent magnet media with linear demagnetizing curve is as follows: where M is the magnetization vector.Substituting (3) in (2) yields The magnetic flux density vector can be represented as the curl of the magnetic vector potential (A): Using ( 4) and ( 5) the following expression is obtained: For each sub-domain, Equation (6) results in Poisson equations for the magnet and slot regions, and Laplace equations for the air-gap and slot-opening sub-domains, as represented below.
In 2-D polar coordinates, the magnetic vector potential and the current density vector have just a component along z, i.e., A = [0, 0, A z (r, θ)] and J = [0, 0, J z (θ, t)].Also, the magnetic flux density vector and the magnetization vector have radial and tangential components as below: Therefore, Equations ( 7)-( 9) are rewritten as

Boundary Conditions
The perpendicular magnetic flux density in two adjacent sub-domains must be equal as mathematically represented as follows: In this equation, B i is the magnetic flux density in sub-domain i, and B i+ is the magnetic flux density in the sub-domain i+.
Also, if there is no current between the two adjacent sub-domains, the tangential components of the magnetic field intensity at the boundary of the two sub-domains are equal; this expression is shown mathematically by Relation (14).
In this equation H i is the magnetic field intensity of the sub-domain i and H i+ is the magnetic field intensity of sub-domain i+.
In both ( 13) and ( 14), n is the perpendicular unit vector to the interface between two adjacent sub-domains.
According to Figure 1, all boundary/interface conditions between sub-domains have been shown from (15) to (25) where α r , β, and δ are respectively the magnet arc per pole pitch ratio, the span angle of slot-openings, and the span angle of slots as shown in Figure 1.

Magnet
Rotor yoke Slot-opening Edges of slot-opening Air-gap Slot-opening Air-gap Slot-opening

Slot
Edge of tooth Slot-opening Edge of tooth

Slot
Stator yoke

Extracting the Fourier Series of the Armature Reaction
To consider the effect of the armature reaction, the current density of each slot should be represented as Fourier series.The current of each phase varies with time and can be represented as Equation (26).
where q is the number of the phases, v shows the order of the harmonics, ω is the angular velocity of the rotor, p is the number of the pole-pairs, γ k = 2π(k − 1)/q is the time offset of the phase kth respect to the first phase.Also I v and θ v are the magnitude and phase offset of vth harmonic, respectively.It is obvious that the relation between the current density in each slot and phase current is dependent on the winding configuration.If the winding configuration is like Figure 2a,b, each slot is considered as one sub-region, like Figure 3a,b.But, if the configuration is similar to Figure 2c, each slot consists of two sub-regions (upper and lower sub-regions), as represented in Figure 3c.
Math.Comput.Appl.2019, 24, 24 7 of 27 In order to complete the Fourier series of each sub-region of a slot, it is necessary to obtain the current density of phases in each sub-region of a slot at a specific time by Equations ( 30) and (31): The current density in each sub-region of a slot can be represented in Fourier series form as in Equation (27).
where J j 0 (t) and J j v (t) are as follows: In order to complete the Fourier series of each sub-region of a slot, it is necessary to obtain the current density of phases in each sub-region of a slot at a specific time by Equations ( 30) and ( 31): For instance, the current density in each sub-domain of a slot could be as Figure 4a,b.
If the winding configuration is as shown in Figure 2a,c, the figure of the current density in each sub-region of slot, at a time instant will be as Figure 4a, and if the winding configuration is as shown in Figure 2b, the figure of the current density in a time instant will be as Figure 4b.Also if the configuration of the winding is similar to Figure 2c, the current density in each sub-region of the slot will be as Figure 4a.In order to complete the Fourier series of each sub-region of a slot, it is necessary to obtain the current density of phases in each sub-region of a slot at a specific time by Equations ( 30) and ( 31): / 2 For instance, the current density in each sub-domain of a slot could be as Figure 4a,b.
If the winding configuration is as shown in Figure 2a or 2c, the figure of the current density in each sub-region of slot, at a time instant will be as Figure 4a, and if the winding configuration is as shown in Figure 2b, the figure of the current density in a time instant will be as Figure 4b.Also if the configuration of the winding is similar to Figure 2c, the current density in each sub-region of the slot will be as Figure 4a.

Extracting the Fourier Series of the Magnetization
In the 2-D polar coordinate system, the magnetization vector only has the radial and tangential components as Equation (32).

Extracting the Fourier Series of the Magnetization
In the 2-D polar coordinate system, the magnetization vector only has the radial and tangential components as Equation (32).
where r and θ are the radial and tangential unit vectors.M r and M θ are the components of magnetization vector which can be represented as Fourier series expansion of Equations ( 33) and (34).
where M rn and M θn respectively are the radial and tangential coefficient of the Fourier series and will be determined according to the magnetization pattern (Table 2).In this paper, only the radial magnetization pattern has been used and represented in Figure 5.
Coefficient of the Tangential Component w−1

Radial Magnetization
where r and θ are the radial and tangential unit vectors.r M and M θ are the components of magnetization vector which can be represented as Fourier series expansion of Equations ( 33) and (34).

Finding the General Solution
The overall format for the general solution in all sub-domains can be represented as Equation ( 35): The general solution not only has the capability to satisfy the related PDE, but must satisfy the boundary conditions of the related sub-domain, especially Equations ( 16), (19), and (24).So the general solutions for sub-domains are as Equations ( 36)- (40).
The general solution of Poisson equation in slot sub-domain will be as follows: ( The particular solution is as follows: The overall format for the general solution in all sub-domains can be represented as Equation (35): The general solution not only has the capability to satisfy the related PDE, but must satisfy the boundary conditions of the related sub-domain, especially Equations ( 16), (19), and (24).So the general solutions for sub-domains are as Equations ( 36)- (40).
The general solution of Poisson equation in slot sub-domain will be as follows: The particular solution is as follows: Also the general solution for the slot-opening sub-domain is The general solution for the air-gap sub-domain is And finally, Poisson equation in the PMs sub-domain has the general solution as Equation (40). where In order to simplify the general solution in PM sub-domain, boundary condition ( 15) has been implemented. where Also, by implementing boundary condition (25), the general solution in slot sub-domain will be simplified as in Equation (45).

Obtaining Integral Coefficients
For implementing boundary condition (17), the correlation technique must be used [31].
From Equation (46), Equation (47) will be deduced. where The solutions of the integrals have been given in the Appendix A. By using correlation technique, Relation (18) must be multiplied by 1 π sin (nθ) and integration on the interval [α − π, α + π], Equation (51) will be obtained.
Equation ( 51) results in Equation ( 52): where The solution of the integral has been given in the Appendix A. Again Equation ( 18) must be multiplied by 1 π cos (nθ), then integration on the interval [α − π, α + π] causes: Equation ( 55) results in Equation (56). 2p−1 where The solution of the integral has been given in the Appendix A.
b so, j 0

Overlapping Winding
In overlapping winding, each slot is divided into two sub-domains as represented in Figure 3c.The upper and lower sub-domains respectively are indicated with slb and slt indices.Hence Equation (69) will be as follows: Also, Equations ( 72) and (75) will be modified as ( 77) and ( 78) respectively.
where R slm = R 2 sl + R 2 so /2 is the radii of middle of the slot which divides it into two equal areas.J j b0 and J j t0 are current densities in the lower and upper sub-domains in a slot.

Flux Density
The air-gap flux density vector is one of the most important quantities required for the calculation of other quantities.For obtaining the air-gap flux density, Equation ( 5) is expanded and Relations (78) and (79) in 2-D polar coordinates are deduced.

Inductances
For calculation of the inductances, just the flux produced by armature is considered.Inductances between phase k and k are obtained as follows: where j = 1, 2, . . ., Q and j = 1, 2, . . ., Q are the indices of the coils.i j is the current of phase k and λ j,j is the flux linked by coil j, which is produced by coil j .If k = k , self-inductance of the phase is calculated.

Back-EMF
In order to calculate the no-load back-EMF of a phase, the permanent magnet flux linked by coils must be calculated by Equation (82).
According to Faraday's law, induced voltage in jth coil obtain by Equation ( 83).
where N t is the number of turns of the coil, ω and α are the angular velocity and rotor position respectively.The total back-EMF of a phase depends on coils connections.

Instantaneous Electromagnetic Torque
Instantaneous torque consists of cogging torque (T cog ), electromagnetic torque (T em ) and reluctance torque (T rel ).
By using Maxwell stress tensor, the instantaneous electromagnetic torque can be obtained as follows: By expanding Equation (85), Relations (86) and (87) will be obtained.
where B a r,PM and B a θ,PM respectively are the radial and tangential flux density components in the air-gap, due to the PMs.Also B a r,AR and B a θ,AR are the radial and tangential magnetic flux density components due to the armature reaction in the air-gap.The parameter R c is the radius of the middle of inner air-gap.

Unbalanced Magnetic Force
The radial and tangential components of the local traction exerted on the rotor surface can be obtained by Maxwell stress tensor as follows: By transforming these local tractions to the Cartesian plane, and summation of the same directions, Equations ( 90)-(93) will be obtained.
Finally, the amplitude of the unbalanced magnetic force can be obtained by Equation (94).

Results
In order to investigate the efficacy of the model, a case study with the parameters listed in Table 3 has been used.Also, the winding configuration for this case study has been shown in Figure 5.

Torque
Instantaneous electromagnetic torque, reluctance torque and cogging torque of the machine have been depicted in Figures 10-12.Both analytic and numeric methods show good agreement which confirms the efficacy of the proposed model.

Torque
Instantaneous electromagnetic torque, reluctance torque and cogging torque of the machine have been depicted in Figures 10-12.Both analytic and numeric methods show good agreement which confirms the efficacy of the proposed model.

Torque
Instantaneous electromagnetic torque, reluctance torque and cogging torque of the machine have been depicted in Figures 10-12.Both analytic and numeric methods show good agreement which confirms the efficacy of the proposed model.

Back-EMF and Inductance
Results of the phase back-EMF and line back-EMF are shown in Figure 13.Again it is shown that both analytic and numeric results have good conformity.
Also, the self and mutual inductances have been depicted in Figure 14.

Back-EMF and Inductance
Results of the phase back-EMF and line back-EMF are shown in Figure 13.Again it is shown that both analytic and numeric results have good conformity.
Also, the self and mutual inductances have been depicted in Figure 14.

Back-EMF and Inductance
Results of the phase back-EMF and line back-EMF are shown in Figure 13.Again it is shown that both analytic and numeric results have good conformity.
Also, the self and mutual inductances have been depicted in Figure 14.

Conclusion
A 2-D analytical magnetic model is presented for brushless synchronous outer rotor machines with surface inset PMs.For this purpose, Maxwell's equations in the form of the Laplace and Poisson equations are solved in predefined sub-domains of the 2-D polar coordinates.The general and particular solutions for each sub-region are presented so that they have the capability to satisfy the governing PDE and related boundary conditions.Finally, by imposing boundary conditions and solving simultaneous linear algebraic equations, all important quantities such as magnetic flux density, electromagnetic torque, UMF, back-EMF, and inductances are calculated and validated by those obtained by FEM.The results of the analytical model show the efficacy of the proposed approach.
Author Contributions: AmirAbbas Vahaj has done formal analysis, investigation, methodology, writing and editing the original draft.Akbar Rahideh has done supervision, conceptualization, project administration and reviewing the original draft.Hossein Moayed-Jahromi has done formal analysis, investigation, methodology, software, and validation.AliReza Ghaffari has done review and editing.

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

Appendix A
To define the left and right side current density in each slot for the two layers non-overlapping concentrated winding topology, we have:

Conclusions
A 2-D analytical magnetic model is presented for brushless synchronous outer rotor machines with surface inset PMs.For this purpose, Maxwell's equations in the form of the Laplace and Poisson equations are solved in predefined sub-domains of the 2-D polar coordinates.The general and particular solutions for each sub-region are presented so that they have the capability to satisfy the governing PDE and related boundary conditions.Finally, by imposing boundary conditions and solving simultaneous linear algebraic equations, all important quantities such as magnetic flux density, electromagnetic torque, UMF, back-EMF, and inductances are calculated and validated by those obtained by FEM.The results of the analytical model show the efficacy of the proposed approach.
(a) According to the geometry and the absence of skewing, the problem is solved in 2-D polar coordinates which means the end effect is neglected.(b) Magnetic vector potential has just axial component which is function of r and θ.Consequently, magnetic flux density has radial and tangential component; i.e., A = [0, 0, A z ], B = [B r , B θ , 0].(c) All materials are isotropic.(d) Rotor and stator back iron have infinite permeability.(e) The edges of the slots and slot-openings have radial direction.(f) The eddy current effect is neglected.

Figure 2 .
Figure 2. Winding topologies (a) single layer alternate teeth wound (b) double layer all teeth wound non-overlapping (c) double-layer overlapping.

Figure 2 .
Figure 2. Winding topologies (a) single layer alternate teeth wound (b) double layer all teeth wound non-overlapping (c) double-layer overlapping.

Figure 3 .
Figure 3. Sub-region division according to the winding configuration: (a) whole slot is considered as one region and belongs to one coil, (b) whole slot is considered as one region and left side and right side of the slot belongs to different coils, (c) whole slot divided into upper and lower sub-regions and each part belongs to one coil.

Figure 3 .
Figure 3. Sub-region division according to the winding configuration: (a) whole slot is considered as one region and belongs to one coil, (b) whole slot is considered as one region and left side and right side of the slot belongs to different coils, (c) whole slot divided into upper and lower sub-regions and each part belongs to one coil.

Figure 3 .
Figure 3. Sub-region division according to the winding configuration: (a) whole slot is considered as one region and belongs to one coil, (b) whole slot is considered as one region and left side and right side of the slot belongs to different coils, (c) whole slot divided into upper and lower sub-regions and each part belongs to one coil.

Figure 4 .
Figure 4. Value of the current density in each sub-region of slot, according to represented winding configurations.

Figure 4 .
Figure 4. Value of the current density in each sub-region of slot, according to represented winding configurations.

Table 3 .
Radial magnetization pattern and its Fourier series components.

7 .
Finding the General Solution in each sub-domain N, U, V, W 100 1

Figure 5 .
Figure 5.The machine topology and winding configuration.

4. 1 .
Flux DensityThe radial and tangential magnetic flux density components due to the open-circuit and armature reaction are respectively depicted in Figures6-9.The numerical results obtained from FEM are shown respectively in each figure and confirm the accuracy of the proposed model.

Figure 5 .
Figure 5.The machine topology and winding configuration.

4. 1 .Figure 6 .
Figure 6.Radial magnetic flux density due to just PMs in the middle of the air-gap, when the rotor position is set to zero.

Figure 7 .Figure 6 .Figure 6 .
Figure 7. Tangential magnetic flux density due to just PMs in the middle of the air-gap, when the rotor position is set to zero.

Figure 7 .Figure 7 .
Figure 7. Tangential magnetic flux density due to just PMs in the middle of the air-gap, when the rotor position is set to zero.

Figure 7 .Figure 8 .
Figure 7. Tangential magnetic flux density due to just PMs in the middle of the air-gap, when the rotor position is set to zero.

Figure 8 .Figure 9 .
Figure 8. Radial magnetic flux density due to just armature winding in the middle of the air-gap, when the current density of phase A is zero, current density of phase B is 4.33 A/mm 2 and phase C is 4.33 A/mm 2 .Math.Comput.Appl.2019, 24, 24 19 of 27

Figure 9 .
Figure 9. Tangential magnetic flux density due to just armature winding in the middle of the air-gap, when the current density of phase A is zero, current density of phase B is 4.33 A/mm 2 and phase C is 4.33 A/mm 2 .

Figure 13 .
Figure 13.Back-electromotive force (EMF) of the first phase and first line.

Figure 13 .
Figure 13.Back-electromotive force (EMF) of the first phase and first line.

Figure 13 .
Figure 13.Back-electromotive force (EMF) of the first phase and first line.

FiFigure 13 .
Figure 13.Back-electromotive force (EMF) of the first phase and first line.

Figure 14 .
Figure 14.Self and mutual inductance of a phase.

4. 4 . 27 Figure 14 .
Figure 14.Self and mutual inductance of a phase.4.4.Unbalanced Magnetic Force (UMF) Unbalanced magnetic forces due to the open-circuit, armature reaction, and both of them have been depicted in Figures 15-17.As evident from these figures, unbalanced forces due to the armature reaction exert considerable forces compared to those of the open circuit.Both analytic and numeric results are shown in Figures 15-17 and confirm the correctness of the proposed model.

Figure 15 .Figure 16 .
Figure 15.Unbalanced magnetic forces just due to the PMs.

Figure 15 . 27 Figure 14 .
Figure 15.Unbalanced magnetic forces just due to the PMs.

4. 4 .
Unbalanced Magnetic Force (UMF) Unbalanced magnetic forces due to the open-circuit, armature reaction, and both of them have been depicted in Figures 15-17.As evident from these figures, unbalanced forces due to the armature reaction exert considerable forces compared to those of the open circuit.Both analytic and numeric results are shown in Figures 15-17 and confirm the correctness of the proposed model.

Figure 15 .Figure 16 .
Figure 15.Unbalanced magnetic forces just due to the PMs.

Figure 16 .
Figure 16.Unbalanced magnetic forces just due to the armature reaction.

Figure 17 .
Figure 17.Total unbalanced magnetic forces exerted on the rotor.

Figure 17 .
Figure 17.Total unbalanced magnetic forces exerted on the rotor.Both analytic and numeric results are shown in Figures 15-17 and confirm the correctness of the proposed model.

Table 1 .
The sub-domains and related symbols.

Table 1 .
The sub-domains and related symbols.

Table 2 .
Radial magnetization pattern and its Fourier series components.

Table 3 .
Parameters used in the case study.