Perspective of Thermal Analysis and Management for Permanent Magnet Machines, with Particular Reference to Hotspot Temperatures

Abstract

Permanent magnet (PM) machines have been extensively used for various applications. Nevertheless, thermal effect, particularly hotspot temperature, not only severely restricts power/torque density but also leads to deteriorations in electromagnetic performance, service life, and reliability. Starting with foundations of PM machines and heat transfer mechanisms, this paper reviews the development of thermal analysis methods over the last thirty years and the state-of-the-art research achievements, and the hotspot temperatures of winding and PM are particularly evaluated. In the overview, various machine losses and cooling techniques are first introduced, which are the essential reasons for temperature rise and the most straightforward way to remove the generated heat. Afterwards, the mainstream thermal analysis techniques, i.e., numerical techniques, lumped-parameter thermal model, and hybrid thermal models, as well as the online electrical parameter-based and thermal model-based temperature monitoring techniques, are reviewed and assessed in depth. In addition, this paper also reviews the analytical thermal modelling methods for winding and PM. Finally, future research trends are highlighted.

1. Introduction

Electrical machines account over 40% of global electricity consumption [1] and play a critical role in both industries and economics. As the evolution of rare-earth permanent magnet (PM) materials [2], the maximum energy products of neodymium iron boron (NdFeB) magnets can achieve towards 500 kJ/m³, Figure 1. Therefore, PM machines can offer high overall efficiency and torque density, and thus, have been extensively used in numerous applications, e.g., aerospace [3], household appliance [4,5], electric/hybrid electric vehicle [6], and wind power generation [7], etc.

Generally, the thermal effect is one of the key factors restricting the increase of power/torque density, especially for the induction and PM machines [8]. On the one hand, the winding insulation will be deteriorated when the local hotspot temperature exceeds the maximum allowed operating temperature. The maximum operating temperature depends on the winding insulation class and is typically <180 °C [9]. Consequently, the degradation of wire insulation could further cause short circuit faults and eventually cause machine failure [10,11]. On the other hand, the thermal deterioration is even more crucial for PM machines since the widely used high energy rare-earth PM material, i.e., NdFeB, is sensitive to overheating. With the increase of PM temperature, the decreasing remanent flux density and the reduction of knee points of PM demagnetization characteristics will reduce output torque capacity and efficiency. More severely, when the local PM hotspot temperature exceeds the maximum allowed operating temperature and the working points operate beyond the knee points of PM demagnetization characteristics, it will cause partial irreversible PM demagnetization. Furthermore, when the PM temperature reaches the Curie temperature, PMs will be fully demagnetized.

Table 1. Typical physical properties of common PM materials [12].

Property	Ferrite	AlNiCo*	SmCo*	NdFeB
$B$ [T]	0.2~0.46	1.1~1.3	0.8~1.2	1.1~1.5
$H_{cJ}$ [kA/m]	210~360	50~150	1300~2400	880~2300
$B{H}_{max}$ [kJ/m3]	6.5~42	35~80	140~260	250~400
${\vartheta }_{max}$ [°C]	180~300	450~525	250	130~230
${\alpha }_{B_r}$ [%/°C]	−0.2	−0.03	−0.06~−0.02	−0.15~−0.1
${\alpha }_{H_{cJ}}$ [%/°C]	0.2~0.5	0.2	−0.4~−0.2	−0.6~−0.4
ρ [Ω⋅cm]×10−6	104	50−80	50–90	110–170

SmCo*: Samarium cobalt. AlNiCo*: Aluminium nickel cobalt.

lists the physical properties of some common PM materials, where $B$ and $H_{cJ}$ are the remanent flux density and the intrinsic coercive force, respectively. $B{H}_{\max }$ is the remanent flux density. ${\vartheta }_{\max }$ is the maximum operating temperature. ${\alpha }_{B_r}$ and ${\alpha }_{H_C}$ are the temperature coefficients of the remanent flux density and the coercive field strength, respectively. $\rho$ is the electrical resistivity.

Although thermal analysis is important for PM machines, it has obtained less attention compared to electromagnetic analysis in the last few decades. This may be attributable to the fact that machine designers may lack a thermodynamic background [13]. More recently, especially in the new century, the importance of the thermal effect has been more well known, and much more effort has been put into the development of thermal analysis [13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33].

Reviewing the evolution of thermal analysis for electrical machines in last few decades, thermal analysis methods, including lumped-parameter thermal models (LPTMs) [34,35,36,37,38,39] and finite element methods (FEMs) [40,41,42,43,44,45], have been employed as early as the 1990s, but were primarily concerned with the induction machines (IM). Meanwhile, the implementation of FEM at that time was complicated and mainly relied on manual programming due to computer hardware limits.

In the first decade of the 21st century, with the popularity of PM machines, the thermal analysis for PM machines garnered much more attention [13,14,15,16,17,18,46,47,48,49,50,51]. The numerical techniques, including electromagnetic-thermal coupled FEM (ETFEM) and computational fluid dynamics (CFD), were maturely programmed in the software and then become one of the most powerful analysis tools [18,47,48,49,50,51,52,53]. LPTMs that were initially modelled for IMs have also been widely applied to PM machines [38,46,47,48,49]. In addition, the LPTM based commercial thermal analysis software Motor-CAD [54] was also developed in 1999 and then further widely employed [13,14,15,16,17,18,19].

From the 2010s, the hybrid thermal models incorporating LPTM with CFD [55,56,57], FEM [58], or the analytical thermal model (ATM) [53,59,60,61,62,63,64,65,66] have been developed to synergize the respective advantages. Meanwhile, several online temperature estimations, including electrical parameter-based [67,68,69,70,71,72,73,74,75,76,77] and thermal model-based techniques [78,79,80,81,82,83,84,85], have also been developed. In addition, machine learning-based neural networks have also been used for thermal monitoring [86,87] recently. In summary, the development of thermal analysis for electrical machines since the 1990s is briefly outlined in Figure 2. In addition, a detailed list of references with remarks is listed in Appendix A.

To date, to help researchers and machine designers be familiar with the thermal analysis of electrical machines, several excellent review papers, e.g., [13,14,15,16,17,18], have been published in the 2000s. However, as aforementioned, numerous new techniques in thermal analysis techniques have been developed over the last 10 years or so which are also worth reviewing and evaluating. More recently, many survey papers have reviewed various cooling techniques specific to different applications, such as wind power generation [19], aerospace [20,21], automotive [22,23,24], advanced manufacturing techniques [28,29], etc., which are summarized in

Table 2. Summary of scopes of recent thermal review papers and this paper.

Ref.	Cooling Methods			Thermal Analysis		Online Temperature Estimation		Remarks
	Passive	Active	Hybrid	Convection	Conduction	Electrical Parameters	Thermal Models
[19]	✓	✓	✓	--	--	--	--	Specific to direct-drive PM generator.
[20]	✓	✓	--	--	--	--	--	Analyses of reliability and cost issues.
[21]	✓	✓	--	--	--	--	--	Review of both electrical machine and power electronics device thermal management.
[22,23,24]	✓	✓	✓	✓	--	--	--	Specific to traction applications and summary of CHT* coefficient equations.
[25]	--	--	--	--	--	✓	✓	Review of electrical parameters, low order and neural network models.
[26]	✓	✓	✓	--	--	--	--	Focusing on hybrid cooling systems on the basis of forced air and liquid cooling.
[27]	✓	✓	--	--	--	--	--	Summary of relationships between different cooling systems and machine power.
[28]	✓	✓	--	--	--	--	--	Focusing on additive manufacturing approaches.
[29]	✓	--	--	--	--	--	--	Review of thermo-conductive materials used for winding insulation system.
[30]	✓	--	--	✓	✓	--	--	Review of uncertain thermal parameters and CHT coefficient equations.
[31,32]	✓	--	--	--	--	--	--	Applications of heat pipes in electrical machine cooling system.
[33]	--	--	--	--	--	✓	--	Online temperature monitoring based on electrical-parameter estimation.
This paper	✓	✓	✓	✓	✓	✓	✓	Comprehensive review of thermal analysis and management with focus on hotspot temperature.

CHT*: Convection heat transfer.

. Obviously, common cooling techniques have been extensively reviewed so far, while the accurate estimation of machine temperatures, particularly hotspot temperatures of winding and PM, are equally or even more important. Nevertheless, the recently developed thermal analysis and temperature estimation techniques, e.g., analytical thermal modelling, the hybrid thermal model, and grey/white thermal model-based online temperature estimation still remain to be reviewed and assessed. Meanwhile, instead of average temperature prediction, hotspot temperature estimation becomes one of the current areas of research, but has not been reviewed so far.

Therefore, with particular reference to hotspot temperatures within electrical machines, this paper reviews the development of thermal analysis and management techniques in the last few decades together with their basic working principles. Meanwhile, mainstream cooling methods, thermal analysis and online temperature estimation methods are also assessed from various perspectives, with specific attention to their effectiveness of hotspot temperature reduction or capabilities of hotspot temperature prediction. The aim is to help readers establish an overall understanding of the state-of-the-art development and challenges, as well as their perspectives in this field.

This paper is organized as follows. Section 2 introduces the foundations of PM machines and heat transfer mechanisms. Section 3 and Section 4 outline the machine loss estimations and cooling methods, respectively. The thermal analysis and online temperature estimation techniques are reviewed and assessed in Section 5 and Section 6, respectively, with particular reference to the hotspot temperatures. Finally, Section 7 concludes the paper and forecasts further research trends.

2. PM Machines and Heat Transfer Mechanisms

PM synchronous machines (PMSMs), as one of the most important brushless electrical machines, have been far more popular in various applications. For the common rotor PM machine, the rotor rotates synchronously with the stator armature reaction field and produces a constant rotor magnetic field by PMs. Thus, the torque can be produced by the interaction of the winding armature reaction and PM fields. The typical radial flux topologies, i.e., surface-mounted PM (SPM) and interior PM (IPM) machines, are illustrated in Figure 3.

Firstly, the SPM machine, as shown in Figure 3a, is the most widely used configuration due to its simple structure, where the PMs are mounted on the rotor surface and adjacent to the airgap. Since the relative permeability of PMs is almost the same as that of the air, there is negligible magnetic saliency and reluctance torque, and only PM torque is produced. For SPM machines, the PM eddy current loss could be significant since the PMs are directly exposed to the armature reaction field. Consequently, large eddy current loss will further increase PM temperature rise and demagnetization risk due to overheating.

Secondly, for the IPM machine, Figure 3b, PMs are embedded in the rotor core, and the leakage flux becomes more significant. Nevertheless, they have a high saliency ratio, i.e., the q-axis inductance is larger than the d-axis inductance, and thereby, they can utilize both PM and reluctance torques to significantly improve the ratio of resultant torque to PM volume. In other words, for the same output torque, the amount of PM materials can be reduced, which is an important advantage of IPM machines when expensive rare-earth magnets are employed. Meanwhile, IPM machines have a better flux-weakening ability and wider speed range since they generally have a higher d-axis inductance compared to that of the equivalent SPM machines. Meanwhile, the design of IPM machines is more flexible since there are more degrees of freedom, which further facilitates their applications. Although the PM loss in the IPM machine is generally lower than that in the SPM machine due to the shielding effect, the PM demagnetization caused by overheating is still a critical issue since the thermal condition on the rotor is usually poor.

In terms of driving currents, the PM machines can be divided into brushless alternating current (BLAC), also called PM synchronous machines (PMSM), and brushless direct current (BLDC) machines, whose back electromotive force (EMF) waveforms are sinusoidal and trapezoidal, respectively. Essentially, BLAC and BLDC machines should be driven under sinusoidal and rectangular phase currents, respectively. However, PM brushless machines may be operated in either BLAC and BLDC mode, although the efficiency, the torque ripple, and the thermal behaviours will be affected [6].

Within PM machines, the thermal energy can be transferred within solids or between the solid and the fluid. There are three types of heat transfer: conduction, convection, and radiation. Figure 4 gives an example of heat transfer paths in a natural cooling IPM machine.

Firstly, conduction heat transfer consists of the propagation of kinetic energy from warm molecules to cool molecules [88]. The conduction heat flux yields Fourier’s law (1), which states that “the heat flux resulting from thermal conduction is proportional to the magnitude of the temperature gradient and opposite to its sign” [89].

$$q_{\mathrm{con}}=-\lambda \nabla T$$

(1)

where $q_{\mathrm{con}}$, $\lambda$, and $\nabla T$ are the conduction heat flux density, the thermal conductivity, and the temperature gradient, respectively.

Secondly, the convection heat transfer is the macroscopic expression of the motion of the particles of the medium and yields Newton’s law of cooling [90], where the convection heat flux density $q_{\mathrm{conv}}$ is defined as:

$$q_{\mathrm{conv}}=h_c\left(T-T_{\infty }\right)$$

(2)

where $h_c$, $T$, and $T_{\infty }$ are the convection heat transfer (CHT) coefficient, the body surface temperature, and the environment fluid temperature, respectively. In PM machines, the convection heat transfers occur within the airgap, the end-cap, between the housing and the environment. Moreover, the widely used cooling techniques, i.e., forced air- and liquid cooling, also mainly rely on enhancing convection heat transfer, which will be introduced in Section 4.

Thirdly, radiation heat transfer is an emission of electromagnetic waves from the objects having above zero absolute temperature. The radiative flux density $q_{\mathrm{rad}}$ can be calculated as:

$$q_{\mathrm{rad}}=\sigma \epsilon T^{4}A$$

(3)

where $\sigma$= 5.669 × 10⁻⁸ W/(m²K⁴) is the Stefan-Boltzmann constant. $\epsilon$ is the emissivity coefficient and depends on the surface material and finish [13,91]. $T$ is the body absolute temperature in Kelvin. In addition, radiation heat transfer usually plays a minor role in electrical machines at normal ambient temperatures and is usually ignored or combined with convection heat transfer.

3. Machine Loss Estimation

During the energy conversion in PM machines, the loss generation is the inevitable parasitic effect, which is also the origin of temperature increases. The main PM machine losses include the stator and rotor iron losses, the copper loss, the PM eddy current loss, and the mechanical loss. The reliable loss estimation is fundamental for accurate thermal estimation. The various losses will be introduced below.

3.1. Stator and Rotor Iron Losses

Iron loss is one of the major losses in the PM machines, and the exploration of accurate iron loss prediction has over one-century history, which was firstly studied by C. Steinmetz [92] in 1892. The Steinmetz model can be written as:

$$P_{\mathrm{Fe}}=k{f}^a{B}^b$$

(4)

where $P_{\mathrm{Fe}}$ is the iron loss, $k$, $a$, and $b$ are Steinmetz coefficients.

With the improved understanding of iron loss, it was found that the generated iron loss can be separated from different origins, i.e., the hysteresis loss $P_{\mathrm{hy}}$, the classical eddy current loss $P_{\mathrm{ec}}$, and the excess loss $P_{\mathrm{exc}}$. To better reflect their different characteristics for engineering convenience, several iron loss models have been proposed and can be divided into two-term models [93,94,95] and three-term models [96,97]. Recently, it has been revealed in [98,99,100,101] that iron loss is also temperature-dependent.

Table 3. Common iron loss models.

Ref.	Iron Loss Model
[94]	$P_{\mathrm{Fe}}=k_{\mathrm{h}}f{B}_m^{\alpha }+k_{\mathrm{ec}}{f}^2{B}_m^2+k_{\mathrm{exc}}{f}^{1.5}{B}_m^{1.5}$
[95]	$P_{\mathrm{Fe}}=k_{h}f{B}_m^{\alpha }+\frac{\sigma d^2}{12T}{{\displaystyle \int }}_0^T{\left(\frac{dB}{dt}\right)}^{2}dt+\frac{k_{\mathrm{exc}}}{T}{{\displaystyle \int }}_0^T{\left|\frac{dB}{dt}\right|}^{1.5}dt$
[96]	$P_{\mathrm{Fe}}=k_{\mathrm{h}}f{B}_m^2+k_{\mathrm{ec}}{f}^2{B}_m^2$
[97]	$P_{\mathrm{Fe}}=k_{\mathrm{h}}\left(f,B_m\right)f{B}_m^2+k_{\mathrm{ec}}\left(f,B_m\right)f^2{B}_m^2$
[98]	$P_{\mathrm{Fe}}=k_{\mathrm{h}}\left(f,B_m,{\vartheta }_{\mathrm{Fe}}\right)f{B}_m^2+k_{\mathrm{ec}}\left(f,B_m,{\vartheta }_{\mathrm{Fe}}\right)f^2{B}_m^2$

lists the common iron loss models, where $k_{\mathrm{h}}$ and $\alpha$ are the hysteresis loss coefficients. $k_{\mathrm{ec}}$ and $k_{\mathrm{exc}}$ are the eddy current and excess loss coefficients, respectively. $B_m$ and $f$ are the amplitude of the flux density and the electrical frequency, respectively. ${\vartheta }_{\mathrm{Fe}}$. is the lamination temperature. In [97], $\sigma$, $d$, and $T$ are the electrical conductivity, the thickness of the lamination, and the time period of the flux density waveform, respectively.

Furthermore, it is reported in [98] that by using the improved thermal-loss coupling method, i.e., considering the temperature dependency of iron loss, the temperature prediction accuracy can be significantly improved. However, it is worth mentioning that the temperature coefficients of hysteresis and eddy current iron losses depend on the flux density and electrical frequency, and thus they need to be determined based on extensive experiments. Moreover, the magnetic properties of silicon steel laminations are also affected by cutting and compression [102,103,104]. It has been demonstrated in the literature that both mechanical and laser cuttings could cause magnetic deterioration, and laser-cut materials have generally worse magnetic properties [102]. Consequently, it gives rise to higher iron loss. Furthermore, the compression of steel lamination caused by shrink fitting also causes a decrease in permeability and an increase in iron loss [104].

3.2. Copper Loss

Copper loss refers to the generated heat due to electrical resistance when the current flows through the wire. For coil-wound electrical machines, the copper loss is a large part of energy loss during the energy conversion, and can be divided into DC and AC components. The DC copper loss is calculated based on the DC resistance, while the AC copper loss is calculated accounting for the skin and proximity effects under the high frequency (HF) AC excitation. The skin and proximity effects represent the influences of the induced eddy currents on the distributions of current density for self-conductor and adjacent conductors, respectively [105]. The analytical calculation of proximity eddy current loss in the conductor can be found in [106,107,108]. Meanwhile, the large flux density variations near the stator slot opening for slotted machines could cause local concentration of AC copper loss for PMSM [109]. Thus, it is reported in [110] that an appropriate arrangement of conductors towards the bottom of the slot can reduce the AC copper loss. Furthermore, it is demonstrated in [88] that the current harmonics caused by the pulse-width modulation (PWM) effect will significantly increase AC copper loss.

Additionally, in comparison to the conventional stranded wire winding, the flat wire (hairpin) winding could increase the slot filling factor and reduce DC copper loss due to its larger cross-section, while the AC copper loss becomes more prominent due to the larger wire diameter [111]. In contrast, the Litz wire could restrain the skin and proximity effects and reduce the AC copper loss [112], while it could increase the DC copper loss due to its drawbacks of low package factor. Therefore, there is a trade-off between the DC and AC copper losses in the windings.

3.3. PM Eddy Current Loss

As is well known, a time-varying magnetic field could induce eddy currents. In an electrically conductive material, the induced eddy current will generate Joule loss. For PM machines during operation, the eddy current loss is generated by asynchronous airgap magnetic fields accounting for the influences of the non-sinusoidal stator winding MMF distribution [113,114], the slotting effect [115,116], the eddy current redistribution field [117,118], and the current harmonics caused by PWM [119,120,121].

Generally, the PM eddy current loss may be relatively smaller compared to the stator iron and copper losses. However, due to the poor thermal condition of the rotor, even a small amount of PM eddy current loss could cause local heat concentration. Meanwhile, as presented in Table 1, the NdFeB PM has a low electrical resistivity and further aggravates its thermal effect and the risk of demagnetization, particularly for the machines having high-speed, high pole number, large slot opening, and/or high power density.

In terms of calculation of PM eddy current loss, the FEM method is a common tool and exhibits a high estimation accuracy. To investigate the physical mechanism of eddy current loss, a huge number of papers regarding analytical calculation can be referred to, such as [116,117,118,119,120,121,122]. In general, the analytical calculation of PM eddy current loss aims to mathematically solve governing equations of magnetic and eddy current fields accounting for the slotting effect, the current harmonics, the eddy current redistribution field, etc. In addition, it is worth mentioning that the distribution of PM eddy current loss is usually non-uniform, which is mainly attributable to the slot opening effect and the armature reaction field. Consequently, the non-uniform loss distribution will give rise to the non-uniform PM temperature increases, and this increases the risk of partial demagnetization due to the hotspot. Detailed analyses of PM temperature distribution will be presented in Section 5.

In terms of reduction of PM eddy current loss, increasing the airgap length is the simplest way, since the high-order abundant field harmonics can be suppressed in the airgap. However, it is not practical for the PM machines requiring high torque density. Moreover, the common PM loss reduction methods include cutting off the complete eddy current paths by using PM axial/circumferential and even radial segmentations [123,124,125,126] or grooving the rotor sleeve [127], suppressing the penetration of MMF sub-harmonics by adding flux barriers on the rotor [128], as well as modifying airgap permeance by adding auxiliary slots from the stator side [129] or shaping the rotor surface [130]. However, there is a trade-off between the manufacturing difficulty and the PM loss reduction.

3.4. Mechancial Loss

Mechanical loss is a general term including windage loss and bearing friction loss. Firstly, windage loss refers to the air friction loss between the rotating rotor and the surrounding air. The windage loss is the energy loss during the acceleration of the surrounding fluid (air) and does not heat the machine but reduces the efficiency [37]. Since the windage loss depends on several factors, i.e., the physical properties of the fluid, the rotor speed, and the geometric parameters of the rotor, it could be more significant for high-speed and/or large volume PM machines.

In general, in electrical machines, the windage losses are generated in the tangential, radial and axial directions. Specifically, the tangential and radial windage losses occur in the airgap and the rotor end-surface due to the rotating rotor. The axial windage loss is generated due to additional forced axial fluid flow, e.g., caused by the axial fan cooling.

The windage losses are usually calculated by the empirical rules, which were developed in the last century and obtained by using the curve fitting technique based on the experiments for Ims. Indeed, they are also applicable for PM machines since the rotors in the PM machines can also be assumed to be a cylinder. The author in [131] summarises the existing empirical rules [132,133,134,135], which will be briefly reviewed below.

The tangential, radial, and axial windage losses $P_{\mathrm{wind},\delta }$, $P_{\mathrm{wind},r}$, $P_{\mathrm{wind},a}$ are calculated as:

$$P_{\mathrm{wind},\delta }=k_1{C}_{\mathrm{f},\delta }{\rho }_{\mathrm{f}}\pi {\omega }^3{r}_o^{4}l$$

(5)

$$P_{\mathrm{wind},r}=0.5C_{\mathrm{f},\mathrm{r}}{\rho }_{\mathrm{f}}{\omega }^3\left(r_o^5-r_i^5\right)$$

(6)

$$P_{\mathrm{wind},a}=k_2{q}_m{\left(\omega r_{\mathrm{o}}\right)}^2$$

(7)

where $k_1$ is the roughness coefficient and equal to 1.0 for smooth surface and 2.5 for axial slotted surface [136]. $C_{\mathrm{f},\delta }$ and $C_{\mathrm{f},\mathrm{r}}$ are the tangential and radial friction coefficients, respectively. ${\rho }_{\mathrm{f}}$ is the fluid density. $\omega$ is the rotor angular velocity. $r_{\mathrm{o}}$ and $r_{\mathrm{i}}$ are the cylinder outer and inner radii, respectively. $l$ is the axial length. $k_2$ is the axial velocity factor and equal to 0.48 for smooth stator and rotor surfaces. $q_m$ is the mass of fluid flow rate within the airgap.

The friction coefficients $C_{\mathrm{f},\delta }$ and $C_{\mathrm{f},\mathrm{r}}$ depend on the cylinder outer radius and the Reynolds number ($Re$). $Re$ is the ratio between the inertia and viscous forces and is used to determine the flow patterns of fluid flow, i.e., the laminar flow at low fluid velocity (small $Re$) or turbulent flow at high fluid velocity (large $Re$). In the tangential and radial directions, the types of fluid flow are determined by the Reynolds numbers $R{e}_{\delta }$ and $R{e}_r$.

$$R{e}_{\delta }=\frac{{\rho }_{\mathrm{f}}\omega r_o\delta }{{\mu }_{\mathrm{f}}}$$

(8)

$$R{e}_r=\frac{{\rho }_{\mathrm{f}}\omega r_o^2}{{\mu }_{\mathrm{f}}}$$

(9)

where $\delta$ and ${\mu }_{\mathrm{f}}$ are the airgap length and the fluid dynamic viscosity, respectively.

Afterwards, the tangential and radial friction coefficients $C_{\mathrm{f},\delta }$ and $C_{\mathrm{f},\mathrm{r}}$ are determined as:

$$\begin{array}{ll}{C}_{f,\delta }=0.515\frac{{\left(\delta /r_o \right)}^{0.3}}{R{e}_{\delta }^{0.5}}\hfill & 500<R{e}_{\delta }<{10}^4\hfill \\ C_{f,\delta }=0.0325\frac{{\left(\delta /r_o \right)}^{0.3}}{R{e}_{\delta }^{0.2}}\hfill & {10}^4<R{e}_{\delta } \hfill \end{array}$$

(10)

$$\left\{\begin{array}{ll}{C}_{\mathrm{f},\mathrm{r}}=\frac{64}{3R{e}_r}\hfill & R{e}_r<30\hfill \\ C_{\mathrm{f},\mathrm{r}}=\frac{3.87}{R{e}_r^{0.5}}\hfill & 30<R{e}_r<3\times {10}^5\hfill \\ C_{\mathrm{f},\mathrm{r}}=\frac{0.146}{R{e}_r^{0.2}}\hfill & \hspace{1em}\hspace{1em}3\times {10}^5<R{e}_r\hfill \end{array}\right.$$

(11)

Secondly, the bearing friction loss $P_{\mathrm{Bearing}}$ is a sum of rolling, sliding, seal, and drag losses, and is mainly dependent on the load and the rotor speed [137].

$$P_{\mathrm{Bearing}}=1.05\times {10}^{-4}{M}_B\times n$$

(12)

$$M_{\mathrm{B}}=\frac{1}{2}{\mu }_{\mathrm{B}}{P}_{\mathrm{L}}{d}_{\mathrm{B}}$$

(13)

where $M_B$ and $n$ are the frictional moment and the rotor speed, respectively. ${\mu }_{\mathrm{B}}$, $P_{\mathrm{L}}$, and $d_{\mathrm{B}}$ are the frictional coefficient, the load on the bearing, and the bore diameter, respectively. The frictional coefficient depends on types of bearing, e.g., ${\mu }_{\mathrm{B}}=0.001$ for a deep groove ball bearing. Meanwhile, the bearing friction loss could also be determined by the no-load deceleration tests and back-to-back drive tests [138,139].

4. Cooling Techniques

The machine cooling technique is a straightforward way to increase the power density of PM machines by removing the generated heat. To date, various cooling techniques [140], including passive cooling, forced air cooling, and forced liquid cooling, have been developed based on the same principles, i.e., enhancement of either conduction or convection heat transfer, as illustrated in Figure 5. Meanwhile, different cooling methods can also be combined to further increase the machine cooling capacitance, which will also be reviewed in this section.

4.1. Passive Cooling

Passive natural cooling is the most common cooling approach for small- and medium-sized electrical machines, where the generated heat dissipates from the housing (frame and end-cap) to the environment via convection heat transfer. For instance, the heat transfer paths within an IPM machine have been shown in Figure 4. Nevertheless, in this case, the high temperature increases of coils and PMs under high power density are difficult to restrain.

Taking account into several practical factors such as cost, size, reliability, and construction complexity, the heat dissipation capability can still be enhanced in a passive way by enlarging the exposed surface of the housing (heat sink/finned housing) [141,142,143], guiding conduction heat transfer [144,145,146], or employing a heat pipe or potting [147,148,149,150].

Firstly, the finned frame in Figure 6 is a useful method to increase the outer surface area and/or the CHT coefficient [141]. It has been demonstrated in [142,143] that the thermal capacity of electrical machines can be improved by optimization of the geometric parameters and number of fins. Although the finned frame/heat sink could decrease the overall temperature levels of the machine, it only has a minor effect on the absolute temperature rises between the housing with the coils/PMs. Thus, its effectiveness in hotspot temperature reduction is limited.

Secondly, the thermal conductors could also be inserted into slots to guide conduction heat transfer from slot to housing (main heat removal path), and thus, restrain hotspot temperatures of active-winding. This is attributed to the higher thermal conductivity of inserted thermal conductors made of aluminium [144] or stator lamination [145] compared to the equivalent thermal conductivity of copper and insulation systems. However, besides the additional induced loss on the metallic conductors, it also brings difficulty in the winding process for a distributed winding configuration. Based on a similar working principle, both global winding [146] and end-winding [147] potting methods can mitigate high-temperature increases in the coils. Since the common potting material, i.e., resin, has a thermal conductivity of 1~3 W/m/K, it is still much higher than the air (0.02~0.03 W/m/K) within the end-cap. It has been demonstrated in the above literature that the hotspot temperature of end-winding can be effectively reduced.

Thirdly, the heat pipe utilizes the phase change of working fluid between the condenser and evaporator to absorb the generated heat. Due to its ultra-high thermal conductivity (>10, 000 W/m/K), a heat pipe has been widely used in different applications [32], including electrical machines [31]. The schematic diagram of various possible applications of heat pipes in electrical machines is shown in Figure 7 [31]. As observed, heat pipes can be assembled on the finned housing, the stator winding, the stator core, the bearing, the rotor, and the shaft. To date, the heat pipes buried in the active-/end-windings are more popular. This demonstrates that the heat pipes inserted in the active-winding [150] and in the end-winding [149] can lead to a 20% reduction of hotspot active-winding temperature and a 20 °C reduction of hotspot end-winding temperature, respectively. However, it needs to be emphasized that the heat dissipation capability of a heat pipe depends on its contact surface area, i.e., higher cooling effectiveness requires a larger contact area.

4.2. Forced Air Cooling

In comparison to passive natural cooling, the forced air cooling approach could increase the CHT coefficients not only around the frame but also within the machine. Normally, the fan is assembled coaxially with the shaft and enhances the fluid flow in the axial direction. The ventilation system can be classified as external fan cooling, internal open and closed fan cooling. Meanwhile, the external and internal (closed) fan cooling could also be employed simultaneously [22]. Figure 8 shows the air flows in different ventilation systems based on an IPM machine.

Firstly, for the external fan cooling in Figure 8a, the fan could accelerate the environmental air flow surrounding the electrical machine and enhance the forced convection heat transfer. The cooling effectiveness depends on the airflow rates, the shape and area of the housing, such as with a finned frame [151]. Secondly, for the internal fan cooling in Figure 8b, the fan could introduce cool air from the external environment into the machine from one side to the other side through filters on the end-plates [152], the airgap, or the stator and rotor ducts. Obviously, since the machine is exposed to the surrounding circumstances, dust, particles, and moisture may pollute electrical machines. Thirdly, to isolate the machine from the surrounding pollutants, the totally enclosed fan cooled (TEFC) ventilation system in Figure 8c is more popular and has been widely used [39,44,141,143,153].

Similarly, the enclosed fan cooling could accelerate the air flows within the machine, and the internal air circulates through stator and/or rotor ducts and the airgap. Furthermore, the circulated air flows could strengthen the heat dissipation from the machine to the environment through the housing. Compared to natural cooling, forced air cooling could reduce overall temperatures within the machine more effectively, while the high temperature rises of winding and PM may not be sufficiently restrained. Meanwhile, the assembled fan could cause additional loss, acoustic noises, and deconstruction risks, especially in high-speed applications.

4.3. Forced Liquid Cooling

In high power applications, e.g., automotive traction machines, the forced liquid cooling technique is more popular due to its effectiveness [18,19]. Generally, compared to the forced air cooling method, the forced liquid cooling system requires additional apparatuses, e.g., high-pressure and scavenge pumps, liquid temperature control units, pressure and flow rate sensors, etc. In addition, the filter may be required when the oil is employed as the coolant. The schematic diagram of a forced liquid system is presented in Figure 9. As observed, the coolant flows from pumps into the electrical machine under high pressure at first. After flowing through the machine, the hotter coolant flows back to the pump. It is also worth mentioning that flow rate and the inlet coolant temperatures play a dominant role in the thermal management of the electrical machine. The flow rate is usually controlled by the throttle and/or bypass valves, while the coolant temperature could be cooled down by a radiator or condenser. In order to achieve satisfactory cooling effectiveness, turbulent flow is desired by increasing the flow rate. Meanwhile, the flow rate can be also adjusted according to the load conditions, especially under overload conditions.

According to whether there is a physical contact between the coolant and the machine components, in other words, whether the coolant needs to be sealed, forced liquid cooling can be divided into two types: the indirect and direct forced liquid cooling approaches. For indirect forced liquid cooling, the coolant does not directly contact the machine components and usually needs to be sealed. Otherwise, it is defined as direct forced liquid cooling in this paper.

4.3.1. Indirect Forced Liquid Cooling

The indirect forced liquid cooling system mainly consists of cooling channels, coolants, and pumps. Firstly, the cooling channels can be grooved or installed within the frame [38,46,48,51,154,155,156], the stator core [157,158,159,160], the slot [161,162], the end-winding [163], and the shaft [164,165,166,167], as presented in Figure 10. Secondly, the liquid coolants could be water, ethylene glycol with water (EGW), or engine oil. Generally, water is preferred due to its high heat capacity. The potential leakage risk is the main issue.

The principle of forced liquid cooling is also based on the enhancement of convection heat transfer. In other words, the cooling capacity is affected by the CHT coefficients and the contact area. The CHT coefficients depend on the geometric parameters and shapes of cooling channels, the physical properties of coolant, and the flow rates [18]. The common frame liquid cooling channel configurations include spiral and axial [54], as shown in Figure 11. It has been reported in [167] that the spiral cooling channel could achieve higher cooling efficiency due to its lower fluid resistance, while the axial channel could achieve a higher contact surface. Similar to the forced air cooling, the frame and stator liquid cooling methods may be insufficient to restrain the hotspot temperature in windings and rotor PMs, which strongly depends on the coolant material, inlet temperature, flow rate, etc.

Furthermore, the cooling channels [161,162] installed in the slots in Figure 10c could significantly reduce winding temperature increases, but the slot filling factor will be sacrificed, and it is not easy to design a compact layout of cooling channels in the slots.

Additionally, the shaft liquid cooling system has been more popular recently, which is especially beneficial with regard to rotor PM temperature reduction [164,165,166], shown in Figure 10d. Meanwhile, it has been demonstrated in [166] that the liquid coolant tends to flow outwardly due to the centrifugal effect during rotation, which could further improve the heat dissipation capability. However, the construction of shaft liquid cooling is more complex, and the mechanical strength of the hollow shaft also needs to be taken into account. In [165], the common direct-through shaft liquid cooling [164,166], shown in Figure 12a, is compared with the recirculating shaft liquid cooling, as shown in Figure 12b. It is also reported in [165] that the recirculating shaft cooling has relatively higher friction loss but a higher CHT coefficient. However, as observed, the construction of the recirculating shaft liquid cooling is more complex.

The cooling channels can also be directly applied to the concentrated end-windings [163], as shown in Figure 10e, while the overhangs of end-winding would be longer than normal to leave the space for channels. It should be emphasised that for the indirect forced liquid cooling techniques, the risks of corrosion and liquid leakage should also be considered.

4.3.2. Direct Forced Liquid Cooling

The direct forced liquid cooling includes immersion cooling [168,169,170] and spray cooling [171,172,173,174,175], as presented in Figure 13. In comparison to indirect cooling, the coolant in a direct cooling system physically contacts specific components, i.e., end-winding and PM, and could provide a more effective temperature reduction due to shorter heat transfer paths.

For the stator immersion cooling in Figure 13a, its cooling principle utilizes high-thermal-conductivity oil to substitute the inner air within the stator part while increasing the CHT coefficient. In this case, the stator and the end-winding are immersed in the oil liquid. The oil could flow through stator ducts from one side to the other side and remove the heat from the stator core and the end-winding via convection heat transfer. However, it needs to be emphasised that the oil sleeve must be specifically designed to avoid oil leakage. In addition, the rotor can also be immersed in the oil liquid, which may lead to larger windage loss during rotor rotation due to the large fluid density of oil [169].

The spray cooling in Figure 13b has become more popular in recent years. An early study in [173] reports that spray cooling with refrigerant coolant provides a better cooling performance compared to immersion cooling. However, the proposed spray cooling is related to the phase change from liquid to gas. Thus, not only convection but also evaporation contributes to heat dissipation. In practice, spray evaporative cooling is still too complicated to be applied in the electrical machine due to the corrosion, pressure, environmental adaptation, compatibility and insulativity of coolants, as well as the additional condenser.

In contrast, spray oil cooling, which is pure convection heat exchange, is favourable due to its high stability and electrical insulation [174], and has also been employed in a Toyota Prius HEV in 2017 [172]. As can be seen from Figure 13b, the spray oil inlets could be located in the frame, the end-cap, or the shaft. The spray cooling is essentially based on the enhanced convection heat transfer between the sprayed cool oil with the end surfaces of end-winding and PM. In [173,175], the influences of oil flow rate, pressure, and the number of nozzles on the cooling efficiency, as well as the CHT coefficient, are investigated. However, it needs to be mentioned that the design of spray oil cooling is complex and exhibits large uncertainty, and the local hotspots in the end-winding may still exist. In general, the higher flow rate nozzles could achieve better cooling performance, while the additional power loss in the pump also needs to be considered. Meanwhile, since the coolant oil directly contacts the electrical insulation system, the long-term corrosion and material compatibility have attracted attention of late. It has been demonstrated in [176] that the different oil-based coolants significantly degrade the mechanical properties of wire coatings, which will also compromise the lifetime of electrical machines. As a summary, the features of different cooling techniques are evaluated in

Table 4. Evaluation of common cooling techniques for electrical machines.

Features	Passive	Forced Air	Forced Liquid
			Indirect			Direct
			Frame or Stator Core	Winding	Shaft	Immersed	Spray oil
Winding hotspot reduction
PM hotspot reduction
Overall effectiveness
Low complexity
Low cost
Maintenance
Corrosion/leakage risks

: Poor. : Medium. : Excellent.

.

4.4. Hybrid Cooling

In some high-power density applications, a single cooling approach may be not enough. Hence, some combined cooling systems have been adopted simultaneously. The hybrid cooling methods can be performed by combining passive and active cooling approaches or by combining two active cooling approaches.

In terms of the combinations of passive and active cooling approaches, the finned frame with forced air cooling is one of the most common methods. In this case, the finned frame can increase not only the contact surface but also the turbulence of fluid, thus increasing the CHT coefficient as well [141,142,143]. Moreover, other passive cooling methods, e.g., additional thermal conductors in the slots, the potting, as well as the heat pipes, are also appropriately employed in forced air or liquid cooling systems [26], which could relatively restrain local overheating. However, it needs to be emphasised that the combined passive cooling components, i.e., inserted thermal conductors, heat pipes, and potting, should be physically connected with the main heat removal path to utilise their temperature difference.

In terms of combinations of the two active cooling approaches, they have been successfully applied to automotive traction machines, which are reviewed in [22,23,24] and are briefly listed below.

• 2012 Tesla Roadster IM: inner forced air + finned housing + outer fan [22].
• 2013 Tesla S60 IM: frame liquid + shaft cooling [24].
• 2014 Porsche Panamera E-hybrid 416: frame liquid + forced air cooling [24].
• GE IPMSM: frame liquid + end-winding spray + rotor cooling [177].
• Zytek PMSM: frame liquid + forced fan cooling [178].

As can be easily found that the combinations of forced liquid and air cooling are widely used in automotive traction machines. However, other possible hybrid cooling methods are also worth further researching accounting for the specific thermal environments and applications. The thermal analysis of electrical machines under different cooling techniques can be implemented based on either LPTM or CFD. The LPTM relies on the CHT coefficients based on the empirical rules [23] and has the advantage of short computation time. However, regarding special and novel cooling systems, CFD and experiments are required to evaluate their cooling effectiveness and to calculate CHT coefficients. The detailed thermal analysis methods will be reviewed in the following section.

5. Thermal Analysis Methods

The thermal effect is one of the most critical factors to be considered at the machine design stage, and the machine designers usually restrict the current density according to the machine size based on previous experience [179,180]. However, it is a very rough indirect thermal analysis and may be not enough to prevent overheating.

Regarding the direct thermal analysis of electrical machines, there are two mainstreams: the numerical techniques, including FEM and CFD, and the conventional high-order LPTMs. Meanwhile, the hybrid thermal models, i.e., FEM + LPTM, CFD + LPTM, and ATM + LPTM, have become more popular in the last decade [181]. The main thermal analysis methods are summarized in Figure 14 and will be reviewed below.

5.1. Numerical Techniques

As is well known, the electromagnetic and thermal phenomena, as well as other most physical problems, could be described by a series of partial differential equations (PDEs), such as the Laplace equation, the Poisson equation, the wave equation, and the Fourier law. The PDEs can be solved by analytical or numerical methods. However, it is not always easy to find analytical solutions, especially for nonlinear problems with complex geometry. In contrast, the FEM becomes one of the most powerful tools for engineers in the fields of electromagnetics, thermodynamics, and fluid dynamics currently. Regarding the working principle of FEM, any large complicated system will be divided into large amounts of small meshed elements. By solving the PDE in every element numerically, a final solution can be obtained in FEM by combining approximated results in all elements.

As early as in the last century, the three-dimensional (3-D) FEM has been used to estimate the thermal behaviours for the turbine generator [182] and IMs [40,41]. In [42], a 3-D thermal FEM is developed to analyze the thermal fields of rotor conductor bars in the IM during starting. Furthermore, the axis-symmetric finite element model was developed in [43] to predict the stator and rotor winding temperatures for IMs under transient conditions. However, due to the limited computer competency, the FEM was mainly programmed manually based on the governing equations.

Since the 21st century, thermal FEM has been successfully included in the commercial electromagnetic analysis packages, including JMAG [183] and ANSYS [184]. Indeed, thermal FEM can accurately calculate conduction heat transfer with arbitrary geometric structures and provide detailed temperature distributions. Meanwhile, it can also be coupled with electromagnetic FEM considering non-uniform loss distributions.

Figure 15 gives examples of steady-state temperature distributions of a 12-slot/10-pole IPMSM calculated by 3-D ETFEM. To date, FEM methods have been extensively used for thermal analysis in PM machines, such as the BLDC machine [51], IPMSMs [46,185,186,187], axial flux PM machines [49,188,189,190], traction PM machines [154,191], and high speed PM machines [136,192,193,194].

On the other hand, the CFD is another numerical technique based on the finite volume method and primarily applied to estimate the fluid fields within and around the electrical machines, as well as the distributions of flow rate and pressure of the coolant in the cooling channels [13]. Before 2000, CFD was used to analyze the fluid field within the end-cap for TEFC IM, accounting for end-winding porosity [44] and geometric changes [45,195]. By using CFD software, the influence of shape of axial fan blade on the heat dissipation for IM is investigated in [196]. In [153], the temperatures for a TEFC IM predicted by the LPTM and the CFD are compared. This demonstrates that the LPTM could achieve high accuracy after thermal parameter tuning. Meanwhile, CFD is also used to estimate CHT coefficients within the end-cap [185,197,198] and the oil-based shaft cooling [164,165,166] for PM machines. It is also employed to estimate the fluid fields for the axial flux PM machines considering the fan cooling [188], the water-cooling [190], and the air-cooling [199] conditions. An example of the CFD simulation result of airflow velocity distribution within a fan-cooling PM machine is presented in Figure 16. In order to simulate the phenomena of turbulent fluid flow, the turbulent models include the one-equation Spalart-Allmaras model, the two-equation k-ε (kinetic energy- dissipation rate of turbulent kinetic energy) and the k–ω (kinetic energy- specific rate of turbulent dissipation) models. Generally, the computational burden increases rapidly with the equation numbers, while the convergence capability is also worse. However, the simple one-equation Spalart-Allmaras model is not appropriate for complex geometry since the turbulent scale is not taken into account in its one-equation governing equation. Therefore, two-equation turbulent models are more popular, where the k-ε model is a widely recognized industrial model and properly used to simulate the turbulent flow away from the wall [22]. In comparison, k–ω model is able to predict the viscous sub-layer and is more suitable for the separation flow perdition.

Overall, the FEM and CFD are useful numerical techniques to analyze conduction and convection heat transfer phenomena with complex geometric structures, respectively. However, they suffer from complicated model setup and extremely high computational burdens. More specifically, the thermal FEM and CFD always require a 3-D model [200,201], which severely increases computational time, especially in transient calculations [60]. Moreover, some components within the PM machines may be not critical in the electromagnetic FEM, e.g., the impregnation, the wire insulation, the slot wedge and liner, the bearing, the housing, etc., but they have a great influence on the thermal behaviours [17]. Consequently, the electromagnetic FEM usually cannot be directly used for thermal analysis. Sometimes, specific knowledge and experience are also required. In contrast, the conventional high-order LPTMs are more preferred due to their easy implementation and short computation time, which will be introduced below.

5.2. Conventional High-Order Lumped-Parameter Thermal Model

In the LPTM, the temperatures, the losses, the thermal capacitances and resistances correspond to the voltages, the currents, the electrical capacitances and resistances in the electrical circuit. The conduction, convection, and radiation heat transfers are represented by the respective thermal resistances, i.e., $R_{\mathrm{cond}}$, $R_{\mathrm{conv}}$, and $R_{\mathrm{rad}}$, which are expressed in (14)–(16), respectively.

In terms of conduction heat transfer within a solid, LPTMs are constituted by either “T-type” [34,60,61,62,63,156] or “I-type” [37,38,65,155] elementary conduction thermal networks, which are used to estimate average ${\vartheta }_{\mathrm{avg}}$ and mid-point ${\vartheta }_{\mathrm{mid}}$ temperatures, respectively. The one-dimensional (1-D) elementary thermal network in the x-direction is shown in Figure 17, where ${\vartheta }_{\mathrm{x},1}$, ${\vartheta }_{\mathrm{x},2}$, and $P$ are the boundary conditions and the power loss, respectively. $R_{\mathrm{x}1}$, $R_{\mathrm{x}2}$, and $R_{\mathrm{xm}}$ with/without superscripts are the conduction thermal resistances. $C$ is the thermal capacitance, as expressed in (17), where $c_p$ and $m$ are the specific heat capacity and the mass of the component, respectively.

$$R_{cond}=\frac{L}{\lambda A}$$

(14)

$$R_{conv}=\frac{1}{h_{c}A}$$

(15)

$$R_{rad}=\frac{1}{h_{R}A}$$

(16)

$$C=c_{p}m$$

(17)

In addition, for multiple-dimensional temperature calculation, the separate thermal networks in all directions are connected to the midpoint, where the heat fluxes in different directions are assumed to be independent. The 2-D or 3-D thermal modelling of each component depends on whether there is heat flow in the third direction. For instance, since there is no heat flow in the circumferential directions for a general cylinder, the 2-D thermal network in the radial-axial plane is appropriate [34]. However, for the thermal modelling of winding, 3-D thermal modelling is necessary since 3-D heat flows exist [202,203].

According to the discretization level, the conventional LPTMs for electrical machines can be classified as:

• Discretized LPTMs [35,36,37,204,205,206].
• Sophisticated LPTMs [46,48,54,59,65,155,207].
• Simplified LPTMs [34,37,38,39,47,154,156,208,209].

Firstly, in the discretized LPTMs, each component within the machine is discretized into several small networks. As an example, Figure 18 shows a 2-D “I-type” based discretized LPTM for a 15 kW TEFC IM [37], which is used to calculate the spatial temperatures for the main components.

Secondly, in the sophisticated LPTM, each component is individually modelled, while the less important component, e.g., rotor core, is lumped into one node to reduce the model complexity. However, for the component of greater thermal significance, e.g., winding, it will be modelled elaborately, e.g., by a 2-D multi-layer model [14,17,65] and a 3-D cylindrical segment [66]. Figure 19 shows a typical sophisticated LPTM built in the Motor-CAD [54] for a 12-slot/10-pole IPMSM [30].

Thirdly, in most cases, simplified LPTMs are more prevalent, where the components having indistinguishable thermal behaviours are merged into one node, such as end-cap and frame. The simplified LPTMs are firstly developed for IMs [34,37,39] and then applied to PMSMs [38,46,47,48,59,60,61,62,63,64,65,66,85,154], and stator PM machines [207,210],. The classical “T-type” and “I-type” based simplified LPTMs for PMSMs are presented in Figure 20.

As can be seen from Figure 18, Figure 19 and Figure 20, all main components within the PMSM need to be involved in any LPTMs to maintain the integrity of the entire thermal network. Meanwhile, the definitions of thermal resistances in (14)–(17) are valid for all LPTMs regardless of the level of discretization. Therefore, there is no substantial difference.

In terms of transient temperature calculation, the thermal network is described by the first-order ordinary differential equation system. In the discrete-time domain, it could be solved directly or approximately by using the well-known Runge-Kutta method [81,82,83,84]. Once the heat fluxes [154] or the temperature rises [155] between two iterative time steps approaches zero, it can be determined that the machine reaches thermal equilibrium.

The characteristics of “I-type” and “T-type” elementary thermal networks are summarized below. The following statements are also applicable to the complete LPTM constituted by the corresponding elementary thermal networks.

Firstly, it is obvious that the “I-type” elementary network has a simple structure to predict the midpoint temperature, while it is developed based on the assumption that the internal loss is generated concentratedly on the midpoint. However, in most cases, the internal loss is generated distributedly throughout the component. Consequently, the “I-type” based LPTM always overestimates the mid-point temperatures due to its concentrated loss input, and the estimation error is proportional to the ratio of the loss in thermal conductivity [211,212]. Thus, this issue is more serious for the active-/end-windings, since they always have low thermal conductivity (<1 W/(m⋅K)) and generate a large copper loss.

Secondly, the “T-type” introduces a negative compensation resistance in each direction to restrain overestimation of 1-D average temperature. However, the estimation performance could be deteriorated significantly in multiple dimensional thermal calculations, especially when the component has a large spatial temperature gradient. This is explained by the fact that in LPTMs, the heat flows are assumed to be independent in multiple dimensions, but the average and hotspot temperature strongly depend on the calculation dimensionality and temperature gradients.

Thirdly, and most importantly, the sophisticated and simplified “T-type” and “I-type” based LPTMs are incapable of estimating the hotspots. The local discretized LPTMs have been applied to assess the slot temperature distributions and hotspots considering the ac copper loss [204] and the back-iron extension [206], as shown in Figure 21. However, the employed “I-type” elementary networks also suffer from inherent defects. Consequently, a considerably large number of nodes (>100) is required to restrain the estimation errors [13,206], but it considerably increases computational burden and modelling difficulty, as shown in Figure 22. Thus, the LPTM loses its superiority of easy implementation.

In summary, sometimes the most concerned maximum temperature cannot be determined or determined accurately in conventional high-order LPTMs. Meanwhile, both the inappropriate calculation of convection coefficients and the misestimation drawbacks would reduce the estimation accuracy, and consequently, give designers the wrong information at the machine design stage. To address this issue, LPTMs have been combined with FEM, CFD, or ATM to utilize the strengths of LPTM and numerical/analytical methods simultaneously, which will be reviewed below.

5.3. Hybrid Thermal Model and Analytical Thermal Modelling

As mentioned above, the hybrid thermal model utilizes the synergies of LPTM with partial CFD or with the FEM/ATM for different purposes. More specifically, regarding convection heat transfer, the partial CFD is used to calculate CHT coefficients for specific regions, which will be further used to determine convection thermal resistances in the LPTM. Meanwhile, regarding conduction heat transfer, the boundary temperatures calculated by LPTM can be employed into FEM or ATM to estimate the temperature distributions for the most important components. The concept and classification of the hybrid thermal models are presented in Figure 23.

The first hybrid thermal model is proposed in [37] and applied to a TEFC IM, where the LPTM and the analytical CFD are combined to calculate the conduction and convection heat transfers (within the end-cap). Furthermore, the partial CFD packages are adopted to estimate only the coolant distribution in the cooling channel and the fluid field within the end-cap for an IM under the steady [41] and transient states [42]. In comparison to the complete CFD modelled for the entire electrical machine, the partial CFD could significantly ease the computational burden [43]. Furthermore, the partial thermal FEM is used to estimate slot temperature distributions and determine the structure of a multi-layer active-winding LPTM for a PMaSynRM [44], where the boundary conditions are calculated by LPTM and then applied to the partial FEM.

In comparison to thermal FEM, the thermal analytical method is more preferable since it can be easily integrated with the LPTM and requires much less computation time with competitive estimation accuracy. Meanwhile, the thermal analytical method is also qualified to estimate temperature distributions and hotspots. This method has garnered much more attention in the last decade [181].

In the hybrid thermal model combining LPTM and ATM, similar to the LPTM+FEM, the boundary conditions are calculated by LPTM, and then they are applied to solve the conduction heat Equation (18).

$$c_i{\rho }_i\frac{\partial {\vartheta }_i}{\partial t}={\lambda }_i\left({\nabla }^2{\vartheta }_i+q_i\right)$$

(18)

where ${\vartheta }_i$ and ${\alpha }_i$ are the temperature and the thermal diffusivity of the solid component (i). ${\nabla }^2$ is the Laplace operator. $q_i$ is the loss density. $c_i$ is the specific heat capacity. ${\rho }_i$ is the mass density. ${\lambda }_i$ is the thermal conductivity. In addition, when the temporal derivative $\partial {\vartheta }_i/\partial t$ equals zero, the conduction heat equation is used to calculate the steady-state temperature distribution.

For active-winding, a simple ATM is proposed in [64] to calculate the active-winding over-temperatures and the feasible current densities at the machine design stage, but the boundary conditions are not considered. In addition, similar to [64] in Figure 24, a circular multi-layered winding model combined with a 2-D active-winding ATM is proposed in [65]. To improve the estimation accuracy, the structure of multi-layered active-winding model is determined by the temperature gradients obtained from the ATM, and need to be “dynamically” adjusted according to different load conditions.

It is worth mentioning that in the multi-layered winding thermal model, either the elliptical multi-layers or the similar triangular multi-layers may suffer from an obvious drawback that the estimation accuracy is sensitive to the thicknesses of virtual impregnation layers due to its extremely low thermal conductivity, i.e., a slight structural change of the virtual impregnation layer will cause a large temperature variation. Consequently, it is challenging to determine an accurate multi-layer structure at the machine design stage due to lack of experimental results.

Both concentrated and distributed end-windings have disparate thermal conditions compared to the active-winding, since end-windings are usually exposed to the inner air within the end-cap. It has been demonstrated in [163] that the temperature of a concentrated end-winding is still ~20 °C higher than that of active-winding under a high-power supply. Therefore, two ATMs for active- and end-windings are modeled separately in [66], where an iterative adjusting process is also developed to overcome the misestimated average temperatures caused by the “T-type” elementary network.

Figure 25 shows an example of steady-state temperature distributions of active- and end-windings for a 12-slot/10-pole IPMSM with concentrated winding [66]. For the active-winding under the copper loss of 9.1 W per slot, the average temperature predicted by the “T-type” network equals 122.9 °C. It is ~5 °C lower than the hotspot temperature of 127.8 °C but is still 8 °C higher than the average temperatures of 115.0 °C calculated by FEM and ATM. The misestimation also occurs in the end-winding. More specifically, a comparison of temperatures predicted by the FEA, the LPTM, and the hybrid thermal model is listed in

Table 5. Comparison of steady-state average and hotspot temperatures of active-/end-windings predicted by main thermal analysis methods [66].

	“T-type” LPTM	Hybrid Thermal Model	Thermal FEA
Active-winding
Average temperature [°C]	122.9	115.0	114.6
Hotspot temperature [°C]	Not predictable	127.6	127.8
End-winding
Average temperature [°C]	123.2	128.2	128.4
Hotspot temperature [°C]	Not predictable	139.7	139.3

, where the detailed information refers to [66]. It can be easily determined that the exemplary results underpin the misestimation issues in conventional LPTM. In comparison, the hybrid thermal model has the similar accuracy compared to thermal FEA.

For PMs, the first hybrid thermal model, LPTM+ATM, is proposed in [59], where a 2-D PM ATM is proposed in the radial-axial plane to compensate for the axial heat flux of PMs. However, the proposed PM ATM cannot determine hotspots since it ignores the temperature distributions in the radial-circumferential plane, where the PM temperature also varies largely. Therefore, the proposed ATM is only validated by the FEM solved under the same governing equation, but may fail to be validated by the ETFEM. Hence, several improved 3-D ATMs have been developed recently, which can predict the steady-state and/or transient spatial PM temperature distributions for both SPMSM with/without rotor sleeve [60,61,62] and IPMSM [63]. Figure 26 shows the steady-state temperature distributions of SPM and IPM calculated by ATMs and 3-D ETFEA [60,63]. In addition,

Table 6. Comparison of steady-state average and hotspot temperatures of SPM and IPMs predicted by main thermal analysis methods [60,63].

	“T-type” LPTM	Hybrid Thermal Model	Thermal FEA
Surface-mounted PM [60]
Average temperature [°C]	164.7	162.8	162.9
Hotspot temperature [°C]	Not predictable	167.7	168.0
Interior PM [63]
Average temperature [°C]	165.5	164.7	165.1
Hotspot temperature [°C]	Not predictable	167.9	168.0

compares the average and hotspot PM temperatures predicted by the FEA, the LPTM, and the hybrid thermal model, where the PM eddy current losses in one PM pole are 7.6 W and 5.1 W, respectively. As observed, the conventional LPTM fails to predict the hotspot temperatures, and the non-uniform distributions of PM eddy current loss also cannot be considered.

Moreover, accounting for the retaining rotor sleeve in high-speed SPM machines, the completed thermal modelling of the rotor is developed in [60], where the different sleeve materials, such as stainless steel or carbon fibre, are taken into account. Based on a 3-slot/2-pole high-speed BLDC machine, the steady-state rotor temperatures of one PM pole in the r-θ plane predicted by the hybrid thermal model and the 3-D ETFEA are shown in Figure 27.

In addition, it is worthy of note that the PM eddy current loss is distributed non-uniformly, which is affected by the slot opening effect, the slot-pole number combinations, and the armature reaction field. Thus, the appropriate mathematical modelling of PM loss is essential. Figure 28 presents the matrix-based post-processing of PM loss and the associated PM temperature estimation method developed in [61].

It can be found that the decrease of PM loss in the radial (r-) direction due to the skin effect is described by a first-order linear equation, while the heat transfer in the circumferential (θ-) direction is ignored due to the periodic symmetry. Furthermore, the standard method of predicting 3-D analytical PM temperature distribution utilizes the multi-slice method based on the axial-radial (z-r) plane [60,63]. The analytical temperature calculation in the 2-D z-r plane can also be degraded to the r-direction, while the axial heat transfer must be considered [62].

As above, there are some common conclusions drawn below.

• The PM temperature distributions and hotspots are strongly affected by the uneven distributions of PM eddy current loss and the boundary conditions.
• The PM end-surface temperature is lower than the PM hotspot due to enhanced convection heat transfer.
• The assumption of uniform PM eddy current loss distribution could cause severe misestimation of PM hotspots.
• For the SPM, the PM loss can be dissipated to airgap via convection. Thus, the SPM hotspots are concentrated on the outer surfaces when the retaining rotor sleeve is not considered.
• For IPM, the PM loss can only be dissipated to the rotor core via conduction. Consequently, the IPM hotspots occur in the upper part inside the PMs due to conduction heat transfer between the PM and the upper rotor core.
• For SPM retained by a rotor sleeve, different sleeve materials have a significant impact on rotor thermal behaviors. Due to the low thermal conductivity of carbon fibre, the hotspot rotor temperature occurs on the adjacent interface between the PM and carbon fibre sleeve. In contrast, the metallic sleeve is beneficial for heat dissipation from the permanent magnet to the airgap. However, the additional eddy current loss induced on the metallic sleeve cannot be ignored.

In summary, Figure 29 summarises the hybrid thermal models and the analytical thermal models specifically for windings, SPM, and IPM, respectively, where the boundary conditions are also highlighted and tabulated in

Table 7. Summary of analytical thermal modelling for windings and PMs.

Component	Ref.	Modelling Plane	Coordinates	Boundary Condition	Analytical Solution
Active-winding	[66,85]	r-θ plane	Polar	Dirichlet (${\vartheta }_{1-4,\mathrm{AW}}$)	Exponential × cosine (sine)
End-winding	[66,85]	x-y plane	Cylindrical	Dirichlet (${\vartheta }_{1-4,\mathrm{EW}}$)	Hyperbolic × cosine (sine)
			Cartesian	Dirichlet (${\vartheta }_{1-4,\mathrm{EW}}$)	Bessel function × cosine (sine)
SPM	[60,61]	r-θ plane	Polar	Robin (${\vartheta }_{\mathrm{SPM},\mathrm{out}}$) Neumann (θ = θ1/2, PM) Dirichlet (${\vartheta }_{\mathrm{SPM},\mathrm{in}}$)	Polynomial + logarithm
		r-z plane	Cartesian	Dirichlet (${\vartheta }_{\mathrm{SPM},\mathrm{out}/\mathrm{in}}$)	Hyperbolic × cosine (sine)
SPM with sleeve	[62]	r-θ plane	Polar	Continuities of temperature and heat flux	Polynomial + logarithm
IPM	[63]	x-z plane	Cartesian	Dirichlet (${\vartheta }_{\mathrm{IPM},\mathrm{out}/\mathrm{in}}$) Neumann (x = w1/2, PM)	Polynomial
		y-z plane	Cartesian	Dirichlet (${\vartheta }_{\mathrm{IPM},\mathrm{out}/\mathrm{in}}$)	Hyperbolic × cosine (sine)

.

5.4. Evalaution and Assessment

Figure 30 and

Table 8. Evaluation of main thermal analysis methods.

Feature	Numerical Techniques	LPTM	Hybrid Thermal Model
Accuracy	High	Relatively low	High
Complex structure	Capable	Medium	Incapable
Computational burden	High	Low	Low
Modelling complexity	High	Low	Low
Estimation of hotspot	Capable	Incapable	Capable
Generality	High	High	Medium

demonstrate the features of the main thermal analysis methods. Indeed, the hybrid thermal model retains the merits of LPTM, i.e., low modelling complexity and computational burden, while also synergizing the strengths of FEM, such as being capable of predicting hotspots with high estimation accuracy. However, it needs to be emphasized that the ATM is mainly appropriate for the components having regular shapes. For example, if the PMs are shaped/grooved and/or have asymmetric geometry considering the minimization of loss [130] and torque ripple [213], numerical techniques are preferred.

5.5. Determination of Uncertain Thermal Parameters

Uncertainties of thermal parameters and properties are the other main challenges in thermal analysis and management. The critical thermal parameters include:

• Material physical properties, including thermal conductivities, specific heat capacities, and mass density.
• Contact thermal resistances of adjacent components caused by imperfect assembling, such as between the frame and the stator core, between the stator lamination and the winding, and between the shaft and the rotor core.
• CHT coefficients within the machine under different cooling conditions.

Firstly, the physical properties of most materials can be found in the datasheet provided by the manufacturer or the data library in the FEA software. The main challenge lies in the modelling of winding, since it is constituted by multiple materials, i.e., the copper, the wire insulation, the impregnation, the slot liner, and/or the slot wedge. In practice, the winding is modelled as a composite solid instead of modelling for all individual wires. The composite solid of winding can be approximated by different models, such as a cylindrical rod [34], a rectangular shape model [37,38], two [39] or multiple layers [13,17,58], a cuboidal model [54], a cylindrical segment [202,203], etc. However, it is challenging to determine the equivalent thermal conductivity. From the perspective of heat transfer, the heat generated in the winding flows in three directions, the radial (r-), the circumferential (θ-), and the axial (z-) directions, as presented in Figure 31.

The radial and circumferential heat must flow through the wire insulation and the impregnation, which usually have extremely low thermal conductivities (<1 W/(m⋅K)). In contrast, the heat transfer along with the axial direction is much easier due to the high thermal conductivity of copper (~400 W/(m⋅K)). Hence, the composite solid has anisotropic thermal conductivities, and the low thermal conductivity in the r- and θ- directions are more critical since they can significantly increase winding temperatures.

Up to now, several approximations have been developed to determine equivalent thermal conductivities. In [34], a simple relationship between the equivalent thermal conductivity and insulation thermal conductivity is proposed. In [37], the Soderberg’s original expression [214] is rearranged and takes account into the effect of airspaces between the slot and the stator lamination. In [14,17], the winding is assumed to be modelled by two components, the inner copper and outer surrounding insulation components. Two regressive equations are developed to estimate the equivalent thermal conductivity of the outer insulation component based on a series of small- and medium-sized induction machines. In addition, the other widely used approximation is proposed by Hashin and Shtrikman (H+S) in [215], which can estimate the equivalent thermal conductivity accounting for two materials (copper and impregnation). More recently, an improved H+S approximation [216] is developed to calculate the composite thermal conductivity with three materials (copper, impregnation, and wire insulation). A summary is given in

Table 9. Summary of formulae of calculating equivalent thermal conductivity of winding.

Ref.	Approximation Formulae
[34]	${\lambda }_{r-\theta }=2.5{\lambda }_{\mathrm{in}}$
[37]	${\lambda }_{r-\theta }=f\left(k_f,{\lambda }_{\mathrm{imp}}\right)$
[17]	${\lambda }_{r-\theta }=0.1076k_f+0.029967$
[14]	${\lambda }_{r-\theta }=0.2425{\left[\left(1-k_f\right)A_{\mathrm{Slot}}{l}_{\mathrm{aw}}\right]}^{-0.4269}$
[215]	${\lambda }_{r-\theta }={\lambda }_{\mathrm{imp}}\frac{\left(1+k_f\right){\lambda }_{\mathrm{cu}}+\left(1-k_f\right){\lambda }_{\mathrm{imp}}}{\left(1-k_f\right){\lambda }_{\mathrm{cu}}+\left(1+k_f\right){\lambda }_{\mathrm{imp}}}$
[216]	${\lambda }_{r-\theta }={\lambda }_{\mathrm{imp}}\frac{\left(1+k_f\right){\lambda }_a+\left(1-k_f\right){\lambda }_{\mathrm{imp}}}{\left(1-k_f\right){\lambda }_a+\left(1+k_f\right){\lambda }_{\mathrm{imp}}}$

, where $k_f$ is the slot filling factor. ${\lambda }_{\mathrm{in}}$, ${\lambda }_{\mathrm{imp}}$, ${\lambda }_{\mathrm{cu}}$, and ${\lambda }_a$ are the thermal conductivities of the insulation, the impregnation, the copper, and the conductor together with the conductor insulation, respectively. $A_{\mathrm{Slot}}$ and $l_{\mathrm{aw}}$ are the slot area and the length of active-winding, respectively.

Secondly, in terms of manufacturing tolerances and surface roughness, the small airspaces caused by the imperfect assembling will significantly increase the temperature rise due to the extremely low thermal conductivity of air (~0.02 W/(m⋅K)). In practice, the airspaces are usually modelled by contact thermal resistances, and the exact values need to be determined based on a dc test by minimizing the predicted and measured temperatures [13,17,37,65].

Ultimately, for natural cooling, the CHT coefficient of the outer housing surface can be approximately equal to 5.95 W/(m²K) [14], which is assumed to be constant for most cases, since the variation of airflow velocity in the environment is usually ignored. For liquid cooling via cooling channels in the housing, the CHT coefficients depend on the shapes and geometric parameters of the channels, the flow rate, and the physical properties of the coolant [14]. Meanwhile, for unconventional cooling systems, e.g., the oil/water-based cooling via hollow shaft [166] and spray oil cooling [173,175], the CHT coefficients need to be determined by CFD software.

More recently, in [77,81], instead of using the empirical rules, the global parameter optimizations are carried out based on a set of comprehensive experiments to determine the thermal resistances and the convection coefficients simultaneously. However, the parameter identifications in [77] cannot reflect the thermal properties within the machine, since the identified LPTM is highly abstract. Meanwhile, the parameter identification in [81] is only concerned with the stator part under standstill with dc supply. Hence, global optimization provides a new perspective of thermal parameter determination and is worth further investigating, while a reliable and robust thermal model is required.

6. Online Temperature Estimation

In comparison to the thermal analysis at the machine design stage, the online temperature estimation during operation is equally or more critical. Normally, PMSMs are designed based on the demanding environment, including the highest ambient temperature, the insufficient cooling capability, and the heaviest load conditions. As a consequence, the machine usually operates with underutilized power since the actual temperatures are unknown.

Against this background, by measuring or estimating the most critical temperatures of windings and PMs, online thermal monitoring can significantly improve the utilization efficiency and prevent insulation failure and PM irreversible demagnetization caused by unexpected overheating. There are usually several thermocouples buried in the active-/end-windings during the winding process. However, the thermocouples increase the number of connections between the machine and the controller, which could reduce the system reliability and increase the complexity [217]. Meanwhile, the inserted thermal sensors may still be incapable of detecting the local hotspots.

On the other hand, it is challenging to assemble thermal sensors in the rotational rotor part during operation. The contact thermal sensors mounted in the rotor part rely on the slip ring, brushes, or hollow shaft for data transmission [218,219,220,221]. The non-contact infrared thermography relies on the visibility of the PM part [58,59,60,61,62,63,72,73,76]. Above all, both contact and non-contact measurement techniques are only appropriate in the laboratory environment and are not suitable in industrial applications since they give rise to additional costs, increase deconstruction risk, and/or are also difficult in terms of maintenance and replacement. Consequently, online temperature estimation methods have been widely investigated in the last decade and can be mainly categorized as the electrical parameter based methods [67,68,69,70,71,72,73,74,75,76,77] and the thermal model based methods [78,79,80,81,82,83,84,85], as shown in Figure 32 and reviewed as follows.

6.1. Electrical Parameter Based Methods

Generally, the electrical parameters, i.e., stator resistance $R_s$ and PM flux linkage ${\psi }_{pm}$, are temperature-dependent (19), and thus, can be used as indicators to reflect average temperatures of winding and PM for online condition monitoring. Since the 2000s, the techniques of online parameter estimation for PM machines have been extensively investigated.

$$R_s\left({\vartheta }_s\right)=R_{s,0}\left[1+\alpha \left({\vartheta }_s-{\vartheta }_0\right)\right]$$

(19a)

$${\psi }_{pm}\left({\vartheta }_{pm}\right)={\psi }_{pm,0}\left[1+\beta \left({\vartheta }_{pm}-{\vartheta }_0\right)\right]$$

(19b)

where $\alpha$ and $\beta$ are the temperature coefficients of copper and PM, respectively. $R_{s,0}$ and ${\psi }_{pm,0}$ are the stator resistance and PM flux linkage at the reference temperature ${\vartheta }_0$. ${\vartheta }_s$ and ${\vartheta }_{pm}$ are the actual coil and PM temperatures, respectively.

In the online electrical parameter estimation, several critical issues need to be emphasized, i.e., the rank-deficient issue and the influence of voltage source inverter (VSI) nonlinearity. The detailed analyses can be found in [33,67]. This paper only gives a brief overview.

Firstly, as observed from the dq-axis steady-state voltage equation in (20), only two parameters can be estimated simultaneously based on one set of dq-axis voltage equations. If the number of unknown parameters is more than the rank of the system, the estimated results could not be converged to the correct values.

$$\left\{\begin{array}{c}{u}_d=R_s\left({\vartheta }_s\right)i_d-{\omega }_r{L}_q{i}_q \hfill \\ u_q=R_s\left({\vartheta }_s\right)i_q+{\omega }_r\left[L_d{i}_d+{\psi }_{pm}\left({\vartheta }_{pm}\right)\right]\end{array}\right.$$

(20)

where $u_d$, $u_q$ are the dq-axis voltages. $i_d$, $i_q$ are the dq-axis currents.${\omega }_r$ is the electrical rotor speed. $L_d$ and $L_q$ are the dq-axis inductances.

On the one hand, some early papers regarding online parameter estimation suffer from the rank-deficiency, i.e., the estimated parameters could not be converged to the correct values [222]. On the other hand, some other papers fix several electrical parameters to estimate others [223], while the estimation accuracy may be deteriorated since both the magnetic saturation and the temperature effect will affect the electrical parameters. These papers will not be discussed.

Secondly, as one of the ill-conditioned issues, the VSI nonlinearity [224,225] severely affects the estimation accuracy of electrical parameters. The state-of-the-art online electrical parameter estimation methods are briefly outlined in Figure 33 and classified by whether magnetic saturation is considered [33].

In the last decade, the electrical parameter based temperature estimation has been extensively investigated. When the dq-axis inductances are known, the stator resistance and PM flux linkage can be estimated separately [69,75] or simultaneously [67,68] to monitor their thermal conditions, while the influences of magnetic saturation and nonlinearity of VSI should be considered. On the other hand, the signal injection methods are also popular, such as current/voltage injection [67,68,69], rotor position offset [70], variable rotor speed [71], and HF signal injection [72,73,74,75]. The schematic diagram of online parameter estimation is shown in Figure 34, where $∆{\omega }_r$ is the variable of the rotor speed. $∆I_d$ and $∆V_d$ are the injected d-axis current and voltage. $∆ \pm {\theta }_p$ are the injected position offsets.

The estimation algorithms include the Adaline neural network, the extended Kalman filter, the model reference adaptive system, the recursive least square, and the particle swarm optimization. It is worth mentioning that a reliable estimation model is the key point. If the estimation model is built appropriately, the different estimation algorithms may mainly affect the convergence speed but play a minor role in the identified values. A review of online parameter estimation and a detailed introduction of estimation algorithms can be found in [33].

Meanwhile, the HF signal injection methods have also been applied to estimate PM temperatures [72,73,74,75]. The HF voltage equations are expressed as:

$$\left[\begin{array}{c}{u}_{d\mathrm{hf}}^r\\ u_{q\mathrm{hf}}^r\end{array}\right]=\left[\begin{array}{cc}{R}_{d\mathrm{hf}}& 0\\ 0& R_{q\mathrm{hf}}\end{array}\right]\left[\begin{array}{c}{i}_{d\mathrm{hf}}^r\\ i_{q\mathrm{hf}}^r\end{array}\right]+\left[\begin{array}{cc}{L}_{d\mathrm{hf}}& 0\\ 0& L_{q\mathrm{hf}}\end{array}\right]\frac{d}{dt}\left[\begin{array}{c}{i}_{d\mathrm{hf}}^r\\ i_{q\mathrm{hf}}^r\end{array}\right]+\left[\begin{array}{cc}0& -{\omega }_r{L}_{q\mathrm{hf}}\\ {\omega }_r{L}_{\mathrm{dhf}}& 0\end{array}\right]\left[\begin{array}{c}{i}_{d\mathrm{hf}}^r\\ i_{q\mathrm{hf}}^r\end{array}\right]$$

(21)

where $u_{d\mathrm{hf}}^r$, $u_{q\mathrm{hf}}^r$, $i_{d\mathrm{hf}}^r$, and $i_{q\mathrm{hf}}^r$ are the dq-axis HF voltages and currents, respectively. $R_{d\mathrm{hf}}$, $R_{q\mathrm{hf}}$, $L_{d\mathrm{hf}}$, and $L_{q\mathrm{hf}}$ are the dq-axis HF resistances and inductances, respectively.

The HF resistance $R_{\mathrm{hf}}$ involves both the stator $R_{\mathrm{shf}}$ and the rotor $R_{\mathrm{rhf}}$ resistances, which are affected by the stator and rotor temperatures ${\vartheta }_s$ and ${\vartheta }_{pm}$, respectively.

$$\begin{array}{ll}{R}_{\mathrm{hf}}& =R_{\mathrm{shf}}\left({\vartheta }_s\right)+R_{\mathrm{rhf}}\left({\vartheta }_{pm}\right)\\ & =R_{\mathrm{shf},0}\left[1+\alpha \left({\vartheta }_s-{\vartheta }_0\right)\right]+R_{\mathrm{rhf},0}\left[1+\beta \left({\vartheta }_{pm}-{\vartheta }_0\right)\right]\end{array}$$

(22)

where $R_{\mathrm{shf},0}$ and $R_{\mathrm{rhf},0}$ are the resistances at the reference temperature ${\vartheta }_0$.

Nevertheless, when $R_{\mathrm{rhf}}$ is used to estimate PM temperatures [72,73], the varying PM magnetization state could cause the change of magnet resistivity due to the magnetoresistive effect, and thus, deteriorate the estimation accuracy. Secondly, $L_{d\mathrm{hf}}$ is another indicator of PM temperature detection [75], since it is associated with the d-axis magnetic saturation and temperature-dependent magnet remanent flux. However, in this case, the information of self- and mutual-magnetic saturations on $L_{d\mathrm{hf}}$ must be known in a prior. In addition, the HF injection-based PM temperature estimations have common drawbacks in that they are inapplicable at the high fundamental frequency, and they also introduce HF losses.

As a summary, the electrical parameter-based temperature estimations are easily integrated into the machine drive system and are capable of tracking the PM ageing effect. Nevertheless, they still suffer from some shortcomings.

Firstly, the temperature coefficients (absolute values) of the copper and the PM are considerably small at 0.39 %/°C and <0.2 %/°C, respectively. When the variations of stator resistance and the PM flux linkage are used to reflect their temperatures (19), the conversion coefficients equal the inverse of the small temperature coefficients. Therefore, the estimated temperatures are very sensitive to measurement errors, VSI nonlinearity, and parameter variations.

Secondly, even if the temperatures could be accurately estimated, they can only represent the average coil and PM temperatures, which may be much lower than the local hotspot temperatures, and thus the online electrical parameter estimation methods are not enough to prevent damages caused by local overheating.

Thirdly, the signal injection-based methods affect the stability of drive systems and could give rise to additional losses and disturbances.

6.2. Thermal Model Based Methods

In comparison, the thermal model-based method is the other mainstream method for online temperature estimation. The adopted thermal model can be classified as white-box and grey-box. In the white-box thermal model, similar to the conventional high-order LPTM, both the topology and the thermal parameters are known [47,63,85]. In contrast, in the grey box, the topology is usually constructed based on experience and knowledge of thermodynamics, while the thermal parameters are unknown and need to be identified [78,79,80,81,82,83,84,85].

Firstly, the thermal model-based temperature estimation should consider the following aspects.

• Ability of predicting transient thermal behaviours.
• Computation time and computing power of the processor.
• Ability to track hotspots (desirable).

6.2.1. White-Box Based Thermal Model

Early in the 2010s, the high-order LPTM based real-time temperature estimation was reported in [47], as shown in Figure 35, where the thermal resistances and capacitances are calculated based on the heat transfer equations and empirical rules. The developed LPTM is represented as state-space equations and solved by using the Euler method for online implementation. However, the high-order LPTM still suffers from the misestimation issue and unknown hotspot, as discussed in Section 5.2.

Thus, more recently, the synergized thermal models [85], combining lumped-parameter and sub-domain thermal models (SDTM), have been developed for a 12-slot/10-pole fractional slot concentrated winding SPMSM in Figure 36, where the detailed machine parameters can be found in [85]. Afterwards, the SDTMs for the active- and end-windings and PM yield respective governing equations (18), and can be solved by utilizing the finite difference method. The interaction between the LPTM and SDTM, as well as the schematic diagrams of SDTMs, are presented in Figure 36c and Figure 37, respectively. The transient temperature distributions and hotspots for these vulnerable components are estimated numerically based on the information on transient boundary conditions. For instance, when the machine operates under 6.8 Nm at 600 r/min and 4.8 Nm at 7200 r/min with the ambient temperature of 20 °C, the corresponding predicted transient temperature distributions of winding and PM at 20 min are shown in Figure 38. It can be seen that the SDTM has a good prediction capability, especially in terms of critical hotspot temperatures.

As a summary,

Table 10. Comparison of characteristics of ATM and SDTM.

Feature	ATM	SDTMs
Governing equation	Conduction heat transfer equation
Boundary condition	Transient/steady-state condition (Table 7)
Loss model	Distributed loss generation
Form of solution	Fourier series with exponential function	Linear equation system
Dependency of accuracy	Accumulated eigenvalues	Mesh grid

compares the characteristics of ATM and SDTM. It can be easily found that their final solutions are the Fourier series multiplying exponential functions and the system of linear equations, respectively. However, there is no essential difference from the perspective of the concept. From this point of view, the analytical thermal models are equally appropriate under online conditions [63], while the analytical solution needs to be updated once the load condition is changed.

6.2.2. Grey-Box Based Thermal Model

However, the white box thermal model, Figure 35 and Figure 36, may be not appropriate in some low-cost applications due to the heavy computation burden caused by their high model order. To address this issue, some low-order grey-box thermal models have been developed subsequently [78,79,80,81,82,83,84] and then became the other popular method for online temperature estimation. In comparison, the low-order LPTMs only consider the most dominant heat transfer paths and the most important components, and thus, they only have two to five nodes to ease the computational burden. An improved low-order grey-box thermal model [84] is presented in Figure 39. It could predict the most critical temperatures of stator core ${\vartheta }_{\mathrm{S}}$, the active-winding ${\vartheta }_{\mathrm{AW}}$, the end-winding ${\vartheta }_{\mathrm{EW}}$, and the PM ${\vartheta }_{\mathrm{PM}}$. ${\vartheta }_{\mathrm{C}}$ and ${\vartheta }_{\mathrm{A}}$ are the reference temperatures the of coolant and ambient, respectively.

In general, the setup procedure of the grey-box thermal model is shown in Figure 40. Firstly, the constructed topology of the grey-box thermal model is essential for estimation accuracy and robustness, i.e., a slight inappropriate structural change would cause a large temperature error. Essentially, the grey-box thermal model should maintain the integrity of the entire heat transfer network of electrical machines to reflect the physical meaning. The inappropriate neglect of components and critical heat transfer paths could cause large estimation errors, especially when the grey-box thermal model is used beyond the training data domain. Furthermore, the redundant thermal resistances in the grey box thermal model could increase the computational burden during the parameter identification. Occasionally, some introduced virtual heat transfer path may negatively affect the model’s robustness.

Secondly, the thermal parameters in the grey-box thermal model are unknown since the heat transfer paths are highly abstract and/or the exact thermal properties or the geometric parameters are unknown. Thus, the thermal resistances are difficult to be calculated analytically in (14)–(17) but are usually identified based on the prior experiments (also known as the training domain). However, since some thermal parameters are varying with the rotor speed, the identification of thermal parameters in the PMSM belongs to a linear parameter varying system. The parameter identification of a linear parameter variation can be implemented by either local approach [78,79,80,81] or global approach [82,83,84].

In the local approach, the identification of a linear parameter varying system is divided into several linear time invariant sub-systems at different constant rotor speeds. For each linear time invariant sub-system, all thermal parameters are kept constant and associated with a constant rotor speed. Ultimately, regarding the entire speed range, all sets of linear time invariant sub-systems are saved into a look-up table. In the global approach, a comprehensive experiment is enough to identify all varying and constant thermal parameters. Thus, it is necessary to construct a mathematical function of varying thermal resistance with rotor speed. However, the highly abstract heat transfer paths would give rise to the difficulty in the construction of functions. Meanwhile, the low-order LPTMs also have drawbacks of the limited degree of freedom due to their abstraction. It has been reported in [84] that this drawback could be mitigated by introducing more degrees of freedom into the thermal resistance definition while not increasing the complexity of the model topology.

Ultimately, in terms of the estimation capability of the grey-box thermal model, after adequate data training, the improved grey-box model (see Figure 39) is employed to predict transient temperatures under two-hour varying load conditions based on a 12-slot/10-pole IPMSM (see Figure 40), as presented in Figure 41. It is obvious that the well-trained grey-box model could achieve an excellent estimation performance, and the maximum temperature errors are kept within 3 °C. However, it should be emphasised that the topology of the grey-box model must be robust, i.e., the abstract grey-box model is able to reflect the physical phenomena to the greatest extent. Meanwhile, the employed identification method also needs to be developed delicately.

6.3. Evaluation and Assessment

In comparison to electrical parameter-based temperature estimation, the thermal model-based methods are robust against measurement errors and parameter variations, and thus have a higher estimation accuracy. However, thermal model based methods require at least one sensor to provide the reference temperature since the temperature increases within the machine are relative values. More importantly, the white-box thermal models, Figure 35, Figure 36, Figure 37 and Figure 38, require accurate information on the geometric parameters and the thermal properties. Even though the exact machine information may be not necessary for the grey-box thermal model in Figure 39, the low-order grey-box thermal models have the drawback that they can only reflect the sensor measured temperatures. If the measured temperatures in the training domain are not the maximum temperatures, the identified model is also incapable of estimating hotspots.

Table 11. Evaluation of main online temperature estimation methods.

Features	Electrical Parameter Based		Thermal Model Based
Requirement of thermal sensors		Unnecessary		At least one sensor
Difficulty of setup or implementation		Easily integrated		Knowledge of thermal properties
Applicability of different PM machines		Salient-/nonsalient-pole		Required individual modelling
Feasibility in entire speed-torque range		Rotor speed dependent		Feasible in the entire range
Capability of multiple temperature estimation		Only coil and PM		All components
Robustness against measurement errors, parameter variations		Sensitive		Robust
Disturbance of drive system		Disturbances (injected signal)		No disturbance
Capability of hotspot tracking		Only average temperatures of coil and PM		Conventional white-box LPTM
				Low-order grey-box (possible)
				Synergized white-box

: Poor. Medium.: Excellent.

summarises the advantages and disadvantages of main online temperature estimation methods.

7. Conclusions and Future Research Trends

With particular attention to the most critical hotspot temperatures, this paper reviews the development of thermal analysis and temperature estimation techniques over the last few decades. The state-of-the-art techniques are also evaluated and discussed. In addition, closely related machine loss estimations and cooling techniques are also outlined.

Currently, with the increasing attention to thermal analysis and thermal management of PM machines, the prospective research trends in this field could be initiated in terms of the following aspects:

1. 

Development of an accurate elementary thermal network.

As previously mentioned, the existing elementary conduction thermal networks more or less suffer from misestimation issues in the multiple-dimensional calculation. Thus, the further development of accurate elementary thermal networks is fundamental but could significantly reduce uncertainties in thermal analysis and temperature estimation.

2. 

Accurate identification of thermal parameters and properties.

The uncertainties of thermal parameters and properties severely affect the robustness and reliability of thermal analysis and temperature estimation. Thus, an accurate method of thermal parameter identification is valuable for both academics and industry.

3. 

Hotspot detection methods for different types of electrical machines.

So far, the recently developed hotspot detection methods are applied to the conventional radial-flux SPM and IPM machines. The specifical thermal modelling and hotspot estimation for other types of electrical machines are also valuable, such as axial flux PM machine, rotor wound synchronous machine, stator PM machine.

4. 

Investigation and application of advanced cooling techniques in different applications of PM machines.

The cooling technique is one of the straightest ways to improve the power density and operational reliability of PM machines. Meanwhile, in different applications accounting for machine configurations, operation conditions, and cost, the most appropriate cooling techniques, including hybrid cooling systems, are also worth investigating.

5. 

Systematic thermal analysis, modelling, and management of the integrated system of PM machine and inverter.

In some compact commercial products, the PM machine and inverter are usually integrated as a unit, and their thermal behaviours interact with each other. Thus, the complete thermal modelling and analysis for the integrated system are also important.