Electromagnetic and Thermal Analysis of Permanent Magnet Synchronous Machines

The increasing interest to permanent magnet (PM) synchronous machines is due to the high torque density and high efficiency that they may exhibit exploiting modern PMs. The three–phase winding is supplied by a current–controlled voltage source inverter, which imposes sinewave currents synchronous with the PM rotor. Such machines are more and more used in several applications, with power rating ranging from fractions of Watts to some Megawatts.


Introduction
The increasing interest to permanent magnet (PM) synchronous machines is due to the high torque density and high efficiency that they may exhibit exploiting modern PMs.The three-phase winding is supplied by a current-controlled voltage source inverter, which imposes sinewave currents synchronous with the PM rotor.Such machines are more and more used in several applications, with power rating ranging from fractions of Watts to some Megawatts.
After a brief introduction on the PM characteristics, this Chapter illustrates the finite element (FE) analysis of the synchronous PM machines.It summarizes the basic concepts dealing with the electromagnetic analysis, and it describes proper analysis strategies in order to predict the PM machine performance.

The permanent magnet machines
In synchronous PM machines, the stator is the same of the induction machines.The rotor can assume different topologies, according to how the PM is placed in it.The machines are distinguished in three classes: surface-mounted PM (SPM) machines, inset PM machine, and interior PM (IPM) machine.Fig. 1(a) shows a cross-section of a four-pole 24-slot SPM machine.There are four PMs mounted with alternate polarity on the surface of the rotor.Fig. 1(b) shows a four-pole inset PM machine, characterized by an iron tooth between each couple of adjacent PMs.Fig. 1(c) shows a four-pole IPM machine, whose rotor is characterized by three flux-barriers per pole.The high number of flux-barriers per pole yield a high rotor anisotropy (Honsinger, 1982).
When the IPM machine is characterized by high anisotropy and moderate PM flux, it is often referred to as PM assisted synchronous reluctance (PMASR) machine.The machine exhibits two torque components: the PM torque and the reluctance torque (Levi, 1984) .

Hard magnetic material (permanent magnet)
The permanent magnets are hard magnetic materials (Bozorth, 1993) .They exhibit a very wide magnetic hysteresis loop, as shown in Fig. 3(a).Once they are magnetized and required to sustain a magnetic field, the PMs operate in quadrant II.Both the intrinsic (dashed line) and normal (solid line) hysteresis loops are drawn in Fig. 3(a), as reported in most PM data sheets.
The intrinsic curve represents the added magnetic flux density that the PM material produces.The normal curve represents the total magnetic flux density which is carried in combination by the air and by the PM (Coey, 1996) .The recoil line of the demagnetization curve is usually approximated by where the residual flux density (or remanence) is B rem and the coercive force is H c (see Fig. 3(a)).The differential relative magnetic permeability of the recoil line is µ rec and is slightly higher than unity.
When the flux density becomes lower than B knee , where the hysteresis curve exhibits a knee in quadrant II, the PM is irreversibly demagnetized and its next operating points move on a locus with lower flux density versus field strength.The minimum flux density in the PMs has to be verified during the magnetic analysis of the machine.Since B knee depends on the temperature, the worst operating condition has to be considered.B knee increases with the temperature, see Fig. 3(b).In order to avoid the irreversible demagnetization of the PM, it is imperative to verify that the minimum flux density in the PM be always higher than the flux density of the knee B knee .
The key properties of some common PM materials are listed in Table 1.Since the operating temperature has a great impact on the PM characteristic, the rate of change of B rem and of H c versus temperature is also reported.

FE analysis pre-processing
In the most of cases, the the magnetic FE analysis is a two-dimensional (2D) analysis 1 .T h e three dimensional effects (e.g., leakage inductance of stator end winding, effect of the rotor skewing, etc.) are not considered here.They have to be computed separately, and added to the 2D field solution.
The simulations are carried out by assuming a planar symmetry: the magnetic field is the same in each section of the machine for a z-axis length equal to the machine stack length L stk .
The current density J and the vector magnetic potential A have only the component normal to the x-y plane, that is, J =(0, 0, J z ) and A =(0, 0, A z ).Therefore the magnetic field strength 1 Let us invite the interested reader to download FEMM (Finite Element Method Magnetics) freeware software by David Meeker.Although FEMM exhibits some limitations (i.e., this means that not all field solution problems can be analysed), it allows the most of electrical machine magnetic field computations to be satisfactorily analyzed.It is very easy to be used, it is organized so as to easily understand how the FE method works.For details about the software, please refer to "FEMM User's Manual"

409
Electromagnetic and Thermal Analysis of Permanent Magnet Synchronous Machines www.intechopen.comH and flux density B vectors have components only in the x-y plane, that is, H =(H x , H y ,0) and B =(B x , B y ,0).Fig. 4(a) shows a cross section of a PM machine in the x-y plane.The relative position between the rotor and the stator reference frames is the mechanical angle ϑ m , which corresponds to the electrical angle ϑ e m = p • ϑ m , where p is the number of pole pairs.

Boundary conditions
In a magnetic field problem, there are four main boundary conditions: Dirichlet: this condition prescribes a given value of the magnetic vector potential A z along a line, typically A z = 0. Thus, the flux lines are tangential to this line.This condition is used to confine the field lines into the domain.
As an example, the Dirichlet boundary condition A z = 0 is assigned to the outer periphery of the stator as shown in Fig. 5(a).
Neumann: this condition imposes the flux density lines to be normal to a line.This condition is a default condition in the field problem.Periodic: this condition is assigned to two lines and imposes that the magnetic vector potential behavior is the same along the two lines, i.e.A z,line1 = A z,line2 .
Anti-periodic: this condition is assigned to two lines and imposes that the magnetic vector potential behavior along one line is opposite to that along the other line, i.e.A z,line1 = −A z,line2 .

Current sources
In addition to the PMs in the rotor, further sources of magnetic field are the stator currents.They are defined by ideal current generators, connected to stator slot surface.In order to impose a prefixed current to each slot, it is convenient to define one generator per each slot Each current is defined by its amplitude and sign.Positive sign means that it is along the direction of z-axis (from the sheet towards reader) and negative sign means that it is opposite to the direction of z-axis (from reader towards the sheet).In magnetostatic analysis, only real part of the current is necessary, while in time harmonic analysis (where symbolic phasor notation is adopted) both real and imaginary components have to be assigned.
Alternatively, a current density can be assigned to the various stator slots.Each slot will be characterized by a different current density, according to the total current flowing in the corresponding slot.The peak value of current density in the slot is defined as where k fill is the slot fill factor and J c is the rms value of the current density in the conductor.
For instance, with k fill = 0.4, the actual current density J c = 6 A/mm 2 corresponds to an equivalent current density Ĵslot = 3.4 A/mm 2 assigned to the stator slot.

Winding definition
The definition of the winding of the electrical machine is extremely important.Fig. 6(a) shows a two-pole three-phase single-layer winding, characterized by two slots per pole per phase.
The arrangement of the coil sides within the slots is described by means of the slot matrix, whose dimension is m × Q, where m is the number of phases (i.e.m=3 in a three-phase machine) and Q is the number of stator slots.The generic element k jq indicates how much the q-th slot is filled by conductors of the j-th phase, where unity means a complete slot fill.
For instance, k jq = 1 means that the q-th slot is completely filled by conductors of the j-th phase; k jq = 0.5 means that only 50% of the q-th slot is filled by conductors of the j-th phase; and k jq = 0 means that no conductor of the j-th phase is in the q-th slot.The element k jq can be either positive or negative sign.The sign refers to the orientation of the coil side.
According to the winding shown in Fig. 6(a), the slot matrix is where the first slot (i.e., q = 1) is the slot in the first quadrant closer to the x-axis, and the other slots follow counterclockwise.
In order to meet the current Kirchhoff law, the element row sum to zero.q 1 2 34 5 6 78 9 1 01 11 2 k a 001 1000 0 -1-10 0 k b -1 -1 0 0 0 0 1 1 0 0 0 0 k c 0 0 0 0 -1 -1 0 0 0 0 1 1 In order to assign the proper current density in each slot, each of them has been defined by a different label, as sketched in Fig. 6(b).The current in the q-th slot can be expressed as where n cs is the number of series conductors per slot, I a , I b , and I c are the currents of the phase a, b, and c, respectively.

Post-processing
After the field solution is achieved, various quantities can be computed from the solved structure.

Flux line plot
The flux lines of the solved structure, shown in Fig. 7(a), give an indication of the flux density that is reached in the various parts of the structure, highlighting the iron saturation.
Similarly, flux density map gives a prompt view of the field solution, allowing to detect gross and striking mistakes in the field problem setting.

Point quantities
In each point of the structure, it is possible to get the value of various electric and magnetic quantities.Fig. 7(b) shows some quantities that can be detected in a generic point P of the structure: magnetic vector potential A z , magnitude and components of flux density vector,

412
Finite Element Analysis -From Biomedical Applications to Industrial Developments www.intechopen.com|B|, B x and B y , magnitude and components of field strength vector, |H|, H x and H y , current density J z .

Magnetic flux
The magnetic flux crossing the air gap is computed by integrating along a line in the middle of the air gap the radial flux density (i.e., the normal component to the line) and multiplying by the z-axis length, that is, The integration line is shown in Fig. 8(a), starting from point S (start) and ending to point E (end).
As an alternative to the line integration, the magnetic flux can be achieved from the magnetic vector potential2 .In this case, it is necessary to compute A z only in the two ends of the line, i.e., in points S and E. The flux results as

Flux linkages
The flux linkage of the phase a is computed referring to the coils arrangement reported in Fig. 6(a).The positive coil sides of the phase a are in the slots 3 and 4, while the negative coil sides are within the slots 9 and 10, as highlighted in Fig. 8(b).Using the slot matrix defined in section 3.4, the flux linkage with the phase a results in where k a,q is the element of the matrix that describes the distribution of the coil sides of the phase a within the q-th slot (see Section 3.4), S slot is the cross-area of the slot.

Magnetic energy
There are three energy quantities that can be computed over the whole structure.The integral of the product of current density by magnetic vector potential Of course, the last integral can be limited in the conductive parts S c of the structure where the current density is not zero.The magnetic energy is HdB dS (8) The magnetic coenergy is BdH dS (9) Even if the magnetic coenergy has not an immediate physical meaning, it is very useful in the magnetic analysis of the machine.Particular care is to define a proper magnetic coenergy density in the PMs which can be defined as

414
Finite Element Analysis -From Biomedical Applications to Industrial Developments www.intechopen.comthat results in a negative value.Alternatively, the magnetic coenergy density in the PMs can be also defined as that is a positive coenergy.Since the lower limit of integration does not affect the computation of the exchange of magnetic coenergy, the "ideal" magnet characteristic is considered and ) is used as the lower limit.Anyway, the rate of change of the magnetic coenergy is the same regardless to the two definitions above, because the difference between the two coenergy densities is the constant quantity 1 2 H * c B rem .
If the linear recoil line ( 1) is adopted for the PM, the magnetic coenergy density results in The second form is particularly advantageous since it could be extended to any material type of the system: hard and soft magnetic material as well as non magnetic material.

Magnetic analysis example
The magnetic finite element (FE) analysis is applied to a synchronous PM machine with Q=24 slots and p=2 pole pairs (Slemon & Straughen, 1980) .Thanks to the machine symmetry, only one pole of the machine is analyzed.

Alignment of the rotor with the stator
The rotor angle is fixed to ϑ e m = 0 when the d-axis (i.e., the rotor reference axis) is coincident with the a-phase axis (i.e., the stator reference axis).This is represented in Fig. 4(b).
Before starting the analysis under load, it is mandatory to adopt the correct reference, that is, to know the correct position of the rotor with respect the stator reference frame.
From the space vector of the PM flux linkage (that is, at no load) the vector angle α e λ is computed, as shown in Fig. 9.If α e λ =0 the rotor is in phase with the stator, otherwise the rotor has to be rotated of a mechanical angle corresponding to the electrical angle of the flux linkage vector, that is:  Then, the d-and q-axis flux linkages are computed adopting the Park transformation.For convention, at no-load there is only d-axis flux linkage, while the q-axis flux linkage is equal to zero.Therefore, at no load, the PM flux linkage is The instant in which the d-axis is aligned to the a-axis represents a particular case.The phase a links the maximum flux due to the PM, and thus Λ m = Λ a .Therefore the PM flux linkage corresponds to the maximum flux linkage of each phase.
The computation can be repeated varying the rotor position.

Computation of the voltage under load
From the flux linkage waveform, it is possible to compute the induced voltage (also called back EMF) waveform, assuming a constant mechanical speed ω m (i.e.electrical speed ω = pω m ).According with the machine convection, it is  The flux linkage waveform is expressed by means of its Fourier series expansion.Then, each harmonic of the series is derived, and all harmonics of EMF are summed together.The final EMF waveform is achieved.An example is reported in Fig. 11.A high distortion in the back EMF waveform is evident in this case.

Flux density in the iron
From the field solution, it is important to verify the maximum flux density in the stator teeth and the maximum flux density in the stator back iron.There are two alternative methods.
The first consists in drawing two circle lines centered in the origin of the axis, one crossing the teeth and the other in the middle of the back iron.The maximum flux density is checked along these two circles.
The second method is based on the hypothesis that the flux density in the stator teeth has essentially radial direction, and the flux density in the stator back iron has essentially azimuthal direction.Therefore, the average flux density in each teeth or in each portion of the back iron can be achieved from the vector magnetic potential in the slots.
As an example, the average flux density in the tooth between slot 1 and 2 is The average flux density in the back iron portion above the slot 1 is since A z = 0 along the external circumference of the stator, where w t is the tooth width and h bi is the back iron height.In order to evaluate the maximum flux density, the computations above have to be repeated for each slots.

Cogging torque
Fig. 12 shows the no-load torque versus the rotor position.This torque is due to the interaction between the rotor PM flux and the stator anisotropy due to the teeth, and it is commonly called cogging torque.In a rotor with identical PM poles equally spaced the number of N p periods during a slot pitch rotation is given by where GCD means Greater Common Divisor.Thus, the mechanical angle corresponding to each period is α τ c = 2π/(N p Q). Higher number of N p periods lower the amplitude of the cogging torque (Bianchi & Bolognani, 2002).Table 2 reports the values of N p for some common combinations of Q and 2p.

Skewing
Skewing rotor PMs, or alternatively stator slots, is a classical method to reduce the cogging torque.In PM machines, the skewing is approximated by placing the PM axially skewed by N s discrete steps (stepped skewing), as illustrated in Fig. 13.When the stepped skewing is adopted, the FE analysis has to be repeated for each of the N s section of the machine, considering the correct angle between the stator and the rotor.
Adopting a stepped skewing, there is also a reduction of the back EMF harmonics and the torque ripple.The analysis is carried out referring to d-q axis component of each electrical and magnetic quantity.

Operation with q-axis stator current
The q-axis current produces a flux in quadrature to the flux due to the PM.When the IPM machine is supplied by q-axis current only, the flux lines go through the rotor and the flux-barrier does not obstruct the q-axis flux, as shown in Fig. 14(a) for a rotor without PM.The q-axis inductance L q exhibits a high value.(Bianchi & Jahns, 2004) .The ratio between the two inductances, i.e. ξ = L q /L d is the saliency ratio (Miller, 1989) .

Magnetic model of the PM synchronous machine
The magnetic model is commonly expressed in the synchronous d-q reference frame.The relationship between the d-and the q-axis flux linkages and currents is non linear, as shown in Fig. 15.It is given by They are single-valued functions, because it is assumed that energy stored in the electromagnetic fields can be described by state functions (White & Woodson, 1959) .Such a model is used for an accurate estimation of the machine performance: to precisely predict the average torque, the torque ripple, the capability to sensorless detect the rotor position.From the d-and q-axis flux linkages, the d-and q-axis voltages are Assuming a linear characteristic of all iron parts, the d-and q-axis flux linkages simplify as Finite Element Analysis -From Biomedical Applications to Industrial Developments www.intechopen.comand the voltage components result in Fig. 16 shows the steady-state 3 vector diagram of the PM synchronous machine in d-q reference frame (Boldea & Nasar, 1999) .

The electromechanical torque
Adopting the FE model of the machine, the torque T can be computed by integrating the Maxwell stress tensor along the rotor periphery: where B g,n and B g,θ are the normal and tangential component of the air gap flux density, and ϑ m is the rotor position (Ida & Bastos, 1992;Jin, 1992;Salon, 1995) .However, assuming the rotor position ϑ m and d-and q-axis currents i d and i q as state variables, the machine torque is also given by where p is the number of pole pairs.W mc is the magnetic coenergy, which is a state function of ϑ m , i d and i q , i.e.W mc = W mc (ϑ m , i d , i q ) (White & Woodson, 1959) .The first term of the second member of ( 23) is labeled as T dq , that is This torque term T dq is slightly affected by the harmonics of the flux linkages and it results to be suitable for the computation of the average torque.The torque ripple is mainly described by the coenergy variation, that is the second term of (23).

Cogging torque
Cogging torque is the ripple torque due to the interaction between the PM flux and the stator teeth.Since the stator currents are zero, it is T dq =0 and, from (23), the cogging torque results to be equal to torque computation (25).A further comparison between predictions and measurements of cogging torque is reported in (Bianchi & Bolognani, 2002) .

Computation under load
Fig. 17 shows the torque behavior versus rotor position of the SPM machine fed by q-axis current only, while d-axis current is zero.Solid line refers to the Maxwell stress tensor computation.The circles refer to the torque computation (23).The dashed line refers to the torque computation T dq , given by (24).As said above, the behaviour of T dq is smooth and close to the average torque.Similar results are found when an IPM machine is considered, as the IPM machine shown in Fig. 1(c).Some torque behaviors and experimental results are reported in (Barcaro et al., 2010A;C) .

Optimizing the torque behaviour
One of the purposes of the FE analysis is to optimize the torque of the machine: high average torque and limited torque ripple (Fratta et al., 1993) .Many strategies exist so as to reduce the 422 Finite Element Analysis -From Biomedical Applications to Industrial Developments www.intechopen.comtorque ripple, and FE analysis is a proper tool to optimize the machine geometry (Bianchi & Bolognani, 2000;Sanada et al., 2004) .
Unusual geometries can be analyzed in order to improve the machine performance.As an example, Fig. 18(a) shows an SPM rotor in which the PMs have been shifted with respect their symmetrical position, so as to reduce the PM flux linkage and the torque harmonics (Bianchi & Bolognani, 2000;2002) .Fig. 18(b) shows the "Machaon" 4 rotor, proposed to reduce the torque ripple in IPM machines (Bianchi et al., 2008;2009) .It is formed by laminations with flux-barriers of different geometry, large and small alternatively under the adjacent poles.

Searching the MTPA trajectory
A proper control strategy is mandatory to achieve high performance of the PM synchronous motors.Fig. 19 shows the key characteristics of the machine as a function of the current vector angle α e i , defined in Fig. 16, for given torque and speed.Normalized parameters are used, that is, unity torque τ pu = 1 and unity speed ω pu = 1 are fixed. 4The name comes from a butterfly with two large and two small wings.

423
Electromagnetic and Thermal Analysis of Permanent Magnet Synchronous Machines www.intechopen.com For a given torque, there is an optimal operating point in which the current is minimum, as shown in Fig. 19(a).Therefore, a maximum torque to current ratio exists.When such a ratio is maximized with respect the current vector angle α e i for any operating condition, the maximum torque-per-Ampere (MTPA) control is achieved (Jahns et al., 1986) .
The flux linkage corresponding to the point of maximum torque is the base flux linkage Λ B .

The flux-weakening operating region
Similarly, the flux-weakening capability of the machine can be investigated starting from the MTPA trajectory, as shown in Fig. 20(b).For each current vector the torque and the d-and q-axis flux linkage components are computed.According to the flux linkage amplitude, Λ, imposing the voltage limit, V n , the maximum electrical speed is computed as ω = V n /Λ.Then, the corresponding mechanical speed ω m = ω/p is related to the electromagnetic torque.
Repeating the computation for different current vector angle α e i , the whole characteristic torque versus speed is achieved.Fig. 21(a) shows the current trajectory, and Fig. 21(b) shows the corresponding torque versus speed curve and power versus speed curve.

Prediction of sensorless capability of PM motors
The technique based on the high-frequency voltage signal injection is used for sensorless rotor position detection of PM synchronous machines at zero and low speed (Ogaswara & Akagi, 1998) .It is strictly bound to the rotor geometry, requiring a synchronous PM machine with anisotropic rotor, e.g. an IPM machine as in Fig. 1(c) or an inset machine as in Fig. 1

(b).
A high-frequency stator (pulsating or rotating) voltage is added to the fundamental voltage, then the corresponding high-frequency stator current is affected by the rotor saliency (Harke et al., 2003;Linke et al., 2003) and information of the rotor position is extracted from current measurement (Consoli et al., 2000;Jang et al., 2003) . 5Actually, the search ends when the torque starts to decrease.
The magnetic model to predict the error signal ε(ϑ e err ) is achieved by a set of finite element simulations carried out so as to compute the d-and q-axis flux linkages as functions of the dand q-axis currents (Bianchi, 2005) .
Then, for a given operating point (defined by the fundamental d-and q-axis currents), a small-signal model is built, defined by the incremental inductances When a high frequency voltage vector is injected along the direction α e v (i.e. the d-and q-axis voltage components are V h cos α e v and V h sin α e v respectively), the small-signal model ( 26) allows to compute the amplitude and the angle of current vector.
The phasor diagram is shown in Fig. 22, including both fundamental components (V 0 and I 0 ) and high frequency components (V h and I h ).Such a study is repeated, varying the voltage vector angle α e v , so as to estimate the rotor position error signal ε.Let I max and I min the maximum and minimum of the high frequency current (computed with the various voltage vector angles α e v ), and α e Imax the angle where I max is found (defined with respect to the d-axis).Then, the rotor position estimation error signal 425 Electromagnetic and Thermal Analysis of Permanent Magnet Synchronous Machines www.intechopen.com is computed as where α e v can be considered as the injection angle and (α e v − α e Imax ) can be considered as the error signal angle ϑ e err .Then, α e Imax corresponds to the angular displacement due to the d-q axis cross-saturation.Fig. 23 compares experimental and predicted results, referring to the inset machine shown in Fig. 1(b), whose rated current is Î = 2.5 A. The pulsating voltage vector technique has been used, adopting a high-frequency voltage with amplitude V h = 50 V, frequency f c = 500 Hz, and a machine speed n = 0 rpm (Bianchi et al., 2007) .The satisfactory match between predictions and measurements, confirms that a PM machine model can be profitably used to predict the sensorless capability of the machine.

Losses in electrical machines
The computation of the losses in an electrical machine is a complex task, sometimes requiring an involved model of the machine even though some uncertainties (e.g. about the material used) prevent a precise estimation.Therefore, in many cases, simpler loss estimations are quite adequate to predict the machine capabilities.

Stator winding losses
Stator resistance is a three-dimensional parameter of the electrical machine.It is generally computed analytically on the basis of the wire diameter and the total length of the winding, including stack length and end-winding length.
The copper conductivity is decreased according to the temperature.Such a temperature is obtained and adjusted after the thermal computation.Sometime an analysis loop is necessary.
Joule losses are computed multiplying the resistance by the square of the rms value of the stator current: P js = 3R s I 2 s .When the current includes time harmonics, the equivalent root-mean-square (rms) current is calculated as: where I 0 is the constant value (if any), and Î1 ... În are the Fourier series coefficients (peak value) of the current waveform.

Stator iron losses
The iron losses of a generic motor consist of the sum of the hysteresis loss, classical eddy current loss and excess loss.Considering a magnetic flux density B varying sinusoidally at the frequency f , the iron loss density is commonly expressed in the following form (Boglietti et al., 2003) as where k hy and k ec are the hysteresis and the classical eddy current constant, and β is the Steinmetz constant, often approximated as β ≃ 2. The k ex is the eddy current excess losses constant.These losses are due to the dynamic losses of the Weiss domains when a variable magnetic field is applied to the magnetic material.The block walls discontinuous movements produce fast Barkhausen jumps and then eddy currents.
These constants should be obtained from material data sheet, but the Epstein frame test does not allow to segregate between the eddy currents due to the classical losses from the eddy currents due to the excess losses.Even if the excess current loss component could be very significant in many lamination materials (Bertotti, 1998) , in (Boglietti et al., 2003) the difficulty to separate the excess losses contribution from the classical eddy current losses has been highlighted.A single eddy current losses coefficient is defined, i.e. an increased k ec , and then applied neglecting the third addendum of (28).
Since the stator iron teeth and the stator back iron operate at different flux density values, the iron loss density has to be computed separately referring to the two parts of the stator.

Tooth iron losses computation with distorted flux density
The equations above hold for sinusoidal flux density variations.When the flux density varies in the iron paths with different waveforms, the computation of the stator iron losses is more complex.
The flux distortions are mainly located in the stator teeth.Let B t the tooth flux density computed as a function of the time (or of the rotor position, when the speed is considered to be constant).An example of the B t waveform is shown in Fig. 24, referring to a four-pole 24-slot machine.Since the machine has two slots per pole per phase, therefore two teeth are considered.It is worth noticing that the waveforms are quite different from sinusoidal waveforms.The behavior of the tooth flux density is expressed by means of Fourier series expansion: where h is the order of the harmonic of the tooth flux density.Only eddy current losses are considered since they are the greatest part of total iron losses due to the flux density harmonics, being proportional to f 2 , see ( 28).
These losses are due to fluctuation of the tooth flux density.Then the tooth eddy current iron loss density is given by: And after some manipulations, it is In (Barcaro et al., 2010B) the effect of the IPM rotor structure on the tooth eddy current iron loss density has been analyzed.

Rotor losses due to MMF harmonics
The discrete location of the coils within the stator slots causes space harmonics of the magneto motive force (MMF) traveling in the air gap.These MMF harmonics move asynchronously with respect to the rotor inducing currents in any conductive rotor parts (Atallah et al., 2000;Shah & Lee, 2006) .The losses in the rotor volume due to the induced currents are given by: where J r is the current density induced in the rotor, A z is the magnetic vector potential, σ is the material conductivity and t is the time.
These losses increase rapidly with the machine size.When the scaling law is apply the flux density B results to be proportionally to the linear quantity l, the curl of B and the induced current density J r to l 2 .Since the volume increases as l 3 , the rotor losses result to be proportionally to l 7 .Although these equations above do not consider skin effects or iron saturation, they highlight how the rotor losses might increase with the size of the machine.The rotor losses phenomenon can be neglected in small size PM machine, it has to be considered compulsorily in large size PM machines.
A practical finite element computation of rotor losses is based on the superposition of the effects: the total rotor losses P rl result as the sum of the rotor losses computed for each harmonic order ν, that is P rl = ∑ ν P rl,ν .To this aim, linear iron is considered assuming an equivalent permeability or freezing the magnetic permeability after a magnetostatic field 428 Finite Element Analysis -From Biomedical Applications to Industrial Developments www.intechopen.comsolution (Bianchi & Bolognani, 1998) .For each ν-th MMF harmonic, the following procedure is adopted: (1) The stator is substituted by an infinitesimal sheet placed at the stator inner diameter D,as shown in Fig. 25.A linear current density K sν (ϑ), sinusoidally distributed in space, is imposed in such a stator sheet: where ν is the harmonic order, ω νr /ν is the speed of such harmonic with respect to the rotor, and Ksν is the peak value of linear current density, achieved from the corresponding MMF harmonic Ûsν .They are: where sgn is equal to +1o r−1 according to whether the harmonic speed is forward or backward the rotor speed (Bianchi & Fornasiero, 2009) .
(2) The circumference is split in a high number of points, e.g.where the phase of the current (i.e.νϑ) is a function of the geometrical position ϑ of the point.
(3) The frequency of the simulation is computed as In each simulation step, the rotor losses due to a single MMF harmonic are computed (Shah & Lee, 2009) .The simulation needs a particular care adopting a two-dimensional analysis.
In each object, electrically insulated by the others, a total current equal to zero is imposed as a further constraint.In addition, when laminations are insulated, a conductivity equal to zero is fixed.The iron conductivity is σ Fe = 3MS/m, and the magnet conductivity is σ PM = 1.16MS/m.Mesh size is chosen according to the penetration thickness.
According to the 12-slot 10-pole PM machine, with a single-layer winding, Fig. 26(a) shows the amplitudes of the MMF harmonics, as percentage of the main harmonic highlighting the presence of a subharmonic of order ν = 1.Fig. 26(b) shows their frequency.Fig. 27 shows the flux lines in the rotor of the 12-slot 10-poles PM machine due to the MMF subharmonic (ν = 1) and other two MMF harmonics of higher order (ν = 7 and ν = 11).

Temperature rises computation
Once the losses of the electrical machine are estimated, it is possible to compute the temperature rise in various parts of the electrical machine.It requires the knowledge of the thermal properties of the material used in the machine, as well as the heat dissipation conditions with the external environment.A two-dimension analysis is considered hereafter.

430
Finite Element Analysis -From Biomedical Applications to Industrial Developments www.intechopen.com

Model of the machine
The FE model of a PM machine is shown in Fig. 28.Each material is characterized by a proper thermal conductivity, and a volume heat generation proportional to its losses.Due to the symmetry of the machine, only a portion of the machine can be analyzed.AB'1")+8, 1)C!"&)61)' Fig. 28.Detail of machine geometry and material properties adopted in the FE analysis.
When the higher losses are located in the stator, a portion formed by one slot and two half-teeth is enough in the analysis.Fig. 28 shows one slot pitch of an interior PM machine.On the stator external surface, there is an Aluminum frame.Between the external surface of the stator lamination and the Aluminum frame, an air foil is added, so as to take into account the roughness and the imperfect contact between the two surfaces.The thickness of such a foil is in the range 0.02-0.05mm.
Within the slot, between the stator lamination and the coil, a slot insulating lining is considered.An insulating separator is added in the middle of these two coil layers to increase the electrical insulation between different phases.

Material thermal properties
The thermal conductivity of the materials are reported in Table 3.As far as the thermal conductivity of the air-gap is concerned, it refers to a fluido-dynamic calculations, using the rotation speed of the machine, as will be described hereafter.
The quality of some materials is dependent on the temperature.Among the others, particular attention is devoted to PMs.Both the residual flux density and the magnetic field of the knee of the B − H curve are reduced, as reported in Fig. 3(b).As said previously the reduction of the magnetic field increases the risk of an irreversible demagnetization of the PMs.

Air gap
The air flow in the air gap is turbulent.The heat transfer is described with an effective thermal conductivity λ gap that is defined as the thermal conductivity that the stationary air should have in order to transfer the same amount of heat as the moving air (Mademlis et al., 2000) .
and ν is the cinematic viscosity of air.Obliviously, if the rotor is at standstill, the air gap conductivity is equal to the stationary value, i.e. λ air .

Slot
The actual slot contains many conductors, insulated each other by means of varnish.However, the drawing of all conductors make no sense.Therefore, the coil winding is modeled as an equivalent homogeneous material characterized by a proper thermal conductivity, as shown in Fig. 29.

Fig. 29. Thermal equivalent model of the slot
The equivalent thermal conductivity is computed as suggested in (Mellor et al., 1991;Schuisky, 1967) .It depends on the ratio between the copper diameter d c and the varnish-insulated diameter d ins .I ti s Finite Element Analysis -From Biomedical Applications to Industrial Developments www.intechopen.comwhere λ ins is the varnish-insulating conductivity, and the F is a multiplier factor (Schuisky, 1967) calculated as F = 37.5 x 2 − 43.75 x + 14 (x = d c /d ins ).Typical value of such a conductivity is 4 to 5 times that of the insulation, for the wire dimensions adopted in practice.

Assigning the loss sources
In each part of the electrical machine in which the losses are generated, a heat generation is imposed as a heat source.Therefore, each region of the model is characterized by specific losses per volume (in W/m 3 ).
The specific losses for the region within the stator slots are computed dividing the Joule losses of the stator winding by the volume of the stator slots.The length L stk of the model has to be considered.
The specific losses for the stator iron region correspond to the stator iron losses divided by the stator iron volume.When the iron tooth losses are quite different from the back iron losses, it is convenient to consider two different regions where to assign the specific losses.In the PM and rotor iron regions, the specific rotor losses are assigned.

Assign the boundary conditions
The temperatures of the winding in the slots, the stator iron and the PMs are referred to the environment temperature, which is fixed to be T env = 0 • C. Therefore, the field solution yields the temperatures rise with respect to the environment temperature.
A natural air convection is considered externally.This yields a thermal convection coefficient equal to 6 W/(m 2 K).However, in the two-dimensional FE model, a higher thermal convection coefficient is set taking into account of the surface increase due to the external length of the frame as respect to the stator stack length, and the presence of fins, which are not detailed in the model of Fig. 28.Typically a factor 3 is used, yielding the thermal convection coefficient to be 18 W/(m 2 K).This is the boundary conditions at the outer surface of the machine frame.
9.5 Thermal computation Fig. 30(b) shows the result of the thermal FE analysis.The steady-state temperature rises are indicated in the slots, in the PMs and in the stator iron.The thermal analysis refers to rated conditions, considering only the Joule losses (standstill operations).
The higher temperature is reached in the slots.The maximum temperature rise reached is about 90 K with respect to the environment temperature.The experimental measurements (Barcaro et al., 2011) confirm the simulated results.According to an external temperature of 20 • C, the winding temperature reaches 110 • C, and the frame temperature 91 • C.

Overload and faulty operating conditions
Overload operations are typically required in many applications, where the operating mode is discontinuous and repetitive accelerations and decelerations are demanded to the electrical machine.Therefore, the electrical machine has to be designed so as to allow temporary overload operations, according to the given cooling system of the system.Particular analysis deals with the faulty operating conditions, in which only a portion of the machine is operating, while the other is disconnected from the power supply.
Hereafter, a dual-three-phase machine is considered in which only half a coils are supplied and the others are open circuited.The machine is fed by a phase current higher than the nominal value, still satisfying the thermal insulating class limit.Referring to the 12-slot 10-pole double-layer PM machine this post-fault strategy has been proposed in (Barcaro et al., 2011) .

Torque under faulty condition
The thermal analysis allows to compute which is the maximum torque that a machine can exhibit according to the limit temperature rise and the fixed cooling system.Some results are presented hereafter referring to the dual-three-phase machine, during faulty operating conditions.Three different coil connections are considered.They are labeled as DL1, DL2, and DL3, according to how the two sets of three-phase windings are placed within the stator.Therefore, according to the thermal analysis, when the dual-three-phase machine is operated under faulty operating conditions, it is possible to reach a torque about 70% of the healthy value overloading the machine and without exceeding the limit temperature.

Impact of the rotor losses
The rotor losses can have a strong impact on the temperature rises of the machine.Fig. 32 shows the temperature map for a PM machine and neglecting and considering rotor losses.
Colors from light blue to dark red show the temperature rise in the machine.In this example, the temperature rises in PMs and Copper increase from 80 K to 100 K, and from 90 K to 98 K, respectively.It is worth noticing how the rotor losses influence the thermal behavior of the two machines.There is an evident increasing of the temperature both in PMs and in the winding, that can cause magnet demagnetization as shown in Fig. 3

Fig. 2
Fig. 2 shows the pictures of a rotor with surface-mounted PMs (a), the lamination of an IPM machine with three flux-barriers per pole (b) and an IPM rotor with two flux-barriers per pole (c).
Fig. 3. (a) Characteristic B-H curve of hard magnetic material, (b) example of Neodymium-Iron-Boron PM including the variation of characteristic with the temperature.

Fig. 4 .
Fig. 4. Electrical machine section in the plane (x, y) and reference frame.

3. 1
Fig. 4(b)  shows the two reference frames.The stator reference frame is characterized by the a-, b-, and c-axis, which refer to the axis of the coils of the three phases.The rotor reference frame is characterized by the d-and q-axis, with d-axis aligned with the PM axis.

Fig. 5 .
Fig. 5. Assignment of boundary condition and current sources.and to assign such a generator to the corresponding slot, as shown in Fig.5(b).Each current is defined by its amplitude and sign.Positive sign means that it is along the direction of z-axis (from the sheet towards reader) and negative sign means that it is opposite to the direction of z-axis (from reader towards the sheet).In magnetostatic analysis, only real part of the current is necessary, while in time harmonic analysis (where symbolic phasor notation is adopted) both real and imaginary components have to be assigned.

411
Fig. 6.Classical representation of the stator winding, and definition of the labels of the stator slots.

Fig. 7 .
Fig. 7. Flux lines and and magnetic quantities detected in the point P of the solved structure.

Fig. 8 .
Fig. 8. Magnetic flux computation from a line integration, and by means of surface integral in the slots containing the coil sides of phase a.

Fig. 10 (Fig. 9 .
Fig. 10(a)  shows one pole of a 24-slot four-pole IPM machine.Fig.10(b)  shows the flux lines at no load, that is, due to the PM only.The flux lines remain in the iron paths of rotor and stator.Some flux lines cross the iron bridges so as to saturate them.

Fig. 10 .
Fig. 10.Flux lines in an IPM machine at no load.

Fig. 11 .
Fig. 11.No load phase and phase-to-phase EMF versus rotor position (at fixed speed).
Fig. 13.Stepped rotor skewing with three modules 5.2 Operation with stator currents Fig. 14.Flux lines in an IPM machine with (a) I q only and (b) I d only.

Fig. 12 Fig. 17 .
Fig.12shows the cogging torque versus rotor position of an SPM machine.Solid line refers to the torque computation by means of the Maxwell stress tensor while the circles refer to the

Fig. 18 .
Fig. 18.Photos of proper solutions to optimize the torque behavior.
Fig. 20.Procedure of changing the current vector angle α e i , searching the MTPA trajectory and defining the flux-weakening capability.

Fig. 22 .
Fig. 22. Phasor diagram with steady-state and high frequency components

Fig. 23 .
Fig. 23.Rotor position error ϑ e err and estimation error signal ε with pulsating voltage injection and inset PM machine (experimental test and prediction).
Analysis of Permanent Magnet Synchronous Machines www.intechopen.com
N p .In each point a prefixed point current I pν is assigned, as shown in Fig.25(b), which is the integral of the distribution of the linear current density over an arc length πD/N p .According to the ν-th harmonic, the maximum current value is computed from the electric loading Ksν ,as Îpν = Ksν πD/N p .The linear current density waveform rotates along the air-gap at the speed ω νr in the rotor reference frame.The points currents are alternating.Using the symbolic notation, in the generic angular position ϑ, the point current is İpν (ϑ)= Îpν e jνϑ (36)

Fig. 30 .
Fig. 30.Machine geometry and simulated temperature rise distribution over environment temperature in healthy mode.

Fig. 31
Fig. 31(d)  shows the temperature map in case of one open-circuited phase according to the DL1 configuration of Fig.31(a).Since the temperature rise of the winding is 52 K (in comparison with 90 K of healthy operating conditions) the operating current could be increased by a factor of √ 90/52 = 1.32.

Fig. 31
Fig. 31(e) shows the temperature map in case of one open-circuited phase according to the DL2 configuration of Fig. 31(b).Fig. 31(f) shows the temperature map in case of one open-circuited phase according to the DL3 configuration of Fig. 31(c).

Fig. 32 .
Fig. 32.Effect of the rotor losses in the machine temperature

Table 2 .
Number N p of cogging torque periods per slot pitch rotation

Table 3 .
Thermal conductivity of the main materialsIts value depends on the relative speed ω m between stator and rotor as well as the air gap width.It is calculated as λ gap = 0.0019 η −2.9084 Re 0.4614 ln(3.3361η)