A Complex Study of Stator Tooth-Coil Winding Thermal Models for PM Synchronous Motors Used in Electric Vehicle Applications

: The operational reliability and high efﬁciency of modern electrical machines depend on the ability to transfer heat in the construction parts of the machine. Therefore, many authors study various thermal models and work on the development of effective heat dissipation. New insights and methods lead to improved techniques for the thermal design of electrical machines. This paper presents an experimentally validated thermal model of a permanent magnet synchronous motor (PMSM) with an improved slot winding model. It also deals with various approaches to homogenization and equivalent material properties of a tooth-coil winding sub-model. First, an algorithm for building a lumped-parameter thermal network (LPTN) of PMSM is described and its properties and problems are discussed. Subsequently, a sub-model of a slot with a winding based on the ﬁnite element method (FEM) is introduced. This sub-model is able to generate different conductor distributions based on probabilistic methods for a speciﬁed ﬁll factor. This allows the veriﬁcation of various homogenization approaches and at the same time it is a tool that automatically calculates thermal resistances for the LPTN.


Introduction
Today, especially with regard to climate challenges, higher energy efficiency and power density of electrical machines are required.This fact brings new challenges in the design of an electrical motor especially for traction applications and electric vehicles (EV).In this day, new types or modified designs of electrical machines are widely analyzed and researched.One of the most utilized types of electrical motors in the automotive industry is a permanent magnet synchronous motor (PMSM).In recent years, research and development of PMSM for automotive electric traction machines have greatly intensified [1][2][3][4].One of the commonly used construction typologies of PMSM in the automotive industry is a machine with tooth-coil winding or a machine with concentrated winding.This brings the advantage of low-cost production, high torque density and high efficiency [5][6][7].On the other hand, the manufacturing process brings new challenges to the thermal design of a PMSM [8].Randomly wound windings mean great uncertainties and the manufactured winding can thus differ significantly from the winding modeled by conventional methods.
The current technical literature divides electric motor thermal analysis and design into two main categories [9]: the thermal design and analyses based on the numerical methods and those based on analytical methods [10,11].The numerical methods include the finite element method (FEM) and computational fluid dynamics (CFD).These methods allow detailed modeling of any machine geometry and the prediction of thermal flow in individual parts of the machine [12][13][14].However, both of these numerical methods are very demanding in the model setup and computational time.Accurate homogenization of the individual motor parts or methods such as model order reduction and cluster computing are then necessary to reduce the computing time [15,16].
An analytical method suitable for the thermal modeling of the electrical machines is the lumped-parameter thermal network (LPTN).The thermal network is based on the machine dimensions and principles of thermal exchange [17,18].One of the unresolved issues of the LPTN approach is the thermal model of the stator winding [16][17][18][19][20].The first problem is that the area of the winding in the slot contains several materials with orders of magnitude different thermal conductivity values.The second problem is that for randomwound windings, it is impossible to determine the exact positions of the conductors in the slot [5].These facts do not allow simple homogenization of the region.Figure 1 shows the differences in the conductor distribution in the slot of one synchronous machine.
Most commonly used methods for modeling of winding slot region are: • Equivalent insulation-This method considers the conductors (copper) located in the middle of the slot region.The insulation and impregnation materials together with the air are modeled as an equivalent insulating material that fills the space around the conductors [16].

•
Layer winding method-This method divides the slot region into several layers of copper and insulation and the model is then able to predict the distribution of temperature in these layers.The temperature of conductors then depends on the positions of the layers [21,22].

•
Cuboidal method-Here, the thermal model of the slot region is based on the cuboidal winding.This method also allows the heat transfer of the winding region into the third (axial) dimension [23].

•
Homogenization method-This method is based on the equivalent thermal conductivity derived from dimensions and material parameters of the slot winding region as a serial connection of the thermal resistances of each material in the slot [10,20].All previously mentioned methods depend on several factors, such as the conductivity of each material, the dimensions of the slot region, the dimensions and the shape of the conductors and the slot fill factor.The fill factor is the ratio between the area of the conductors and the slot area.A higher fill factor means a higher power density of the electrical machine and a better thermal dissipation from the slot area.For various winding types, the slot fill factor ranges from 40% to 85%.The simplest winding types are characterized by round conductors, a random-wound winding process with great uncertainties, low manufacturing requirements and therefore low cost.These types are suitable for mass production, but they have a lower fill factor around 40% [24].This is the main reason why they are used for the high-power traction machines with lower fill factors.These machines require high efficiency, high power density and low-cost manufacturing [25][26][27].Another tooth-coil winding type with round conductors is compressed winding.Thanks to the high packing density of conductors, this winding type achieves a higher fill factor approximately 65% [24].The winding with the highest fill factor is layer winding, where the conductors are uniformly arranged in layers [28].
The main goal of this paper is to investigate the homogenization and estimation of the equivalent material properties of the sub-model of the tooth-coil winding.As mentioned above this is currently one of the main topics discussed in the thermal modeling of an electrical machines.Our approach uniquely combines LPTN and FE models, achieving an error of less than 5% compared to the measurement.The LPTN model is fast, modular and allows easy redefinition and addition of individual parts for new machine topologies.In order to increase the accuracy of the winding temperature calculation, the FE model of the slot with winding is presented.This model is based on a random conductor generator that simulates the actual process of manufacturing of randomly wound windings.The FE model calculates the thermal resistances for the LPTN model and further allows a more accurate calculation of temperature distribution after the boundary conditions are known from the LPTN model.Our hybrid model is then verified by measuring temperatures on a traction high-speed synchronous motor.
An important parameter that is mentioned throughout the text is the equivalent thermal conductivity k eqv .This parameter expresses in one number the thermal conductivity of a certain part of the machine.There are several possible methods for calculating this conductivity and they are presented in the following sections.
The work presented in this paper is organized into six chapters.A flowchart of the research is shown in Figure 2 and it also corresponds to the paper's structure.Each chapter is focused on the investigation and description of the individual part of a main topic.Section 1 describes the state-of-the-art thermal modeling of electrical machines; this section also defines the main research problem in thermal modeling and presents the goal of this paper.Section 2 describes the software algorithm of the modular thermal model of the PMSM.Section 3 presents the algorithm of the finite element methods with a random conductor distribution generator and also shows the relation between slot fill factor and the range of thermal resistances of the winding.Section 4 shows experimental validation of calculated results for various methods for estimating equivalent thermal resistances presented in Sections 2 and 3. Discussion of results obtained by measurement and presented thermal models is presented in Section 5. Investigated results are summarized in Section 6.

Software Algorithmization of Thermal Modeling
When designing a machine, it is necessary to quickly assess various thermal designs, and LPTN models can serve well for this.To be as efficient as possible, the model should automatically generate a thermal network based on the input parameters.Furthermore, it is necessary to be able to implement new machine topologies and new research results into the model as easily as possible.This is achieved through the concept of sub-models, which divide the electrical motor into individual layers from shaft to housing.These layers are shown in Figure 3.This model also allows extended, more accurate thermal calculations of individual sub-models after the global model is calculated and the temperatures of the surrounding components are known.The software is built to be as clear as possible for the user and there is also a possibility to connect the results of electromagnetic design procedure to the input file.The software is currently limited to steady-state analysis of a PMSM with standard rotor and radial magnetic flux.The calculation procedure is shown in Figure 4.  First, all necessary data are set or loaded, such as the lengths of discretization elements, dimensions of machine parts, ambient temperatures, etc.In the next step the number of a nodes per component is set, based on the length of the elements and dimensions of the component.Next the global thermal network is defined, which requires an estimation of the convection heat transfer coefficients.When this is done, the sub-models are defined.The sub-model creation is shown in Figure 4.For each homogeneous layer, thermal resistances are calculated from input parameters and thermal connections are established.Finally, a global matrix is assembled, the vector of heat losses for individual nodes is set and steady-state analysis can be computed.If necessary, further calculations are performed and all resulting temperatures are processed into a result table.

Lumped Thermal Models of Winding
As already mentioned, the winding is one of the parts of the machine that is very problematic from the point of view of thermal modeling.Due to uncertain distribution of conductors in the slot and imperfect impregnation, which leaves air cavities with poor thermal conductivity, it is very difficult to correctly predict the thermal distribution of the winding by analytical methods.In our software the general winding is modeled as the element shown in Figure 5, where R AX are thermal resistances in the axial direction, R RAD are resistances in the radial direction towards the yoke and towards the air gap and R TAN is the tangential resistance towards the stator tooth.In the axial direction, only the area of copper, which can be easily analytically described, is used to determine the resistance.However, it is difficult to determine the resistances in the radial and tangential direction, for which separate models are required for each type of winding.For windings with rectangular conductors the analytical model is quite accurate because positions of conductors and layers of insulation are known.For this reason, it is possible to create a detailed thermal network with a node in each conductor.From its solution the maximum temperature and the heat flow in each direction can be determined.The temperature difference and heat flux are then used to calculate the thermal resistances shown in Figure 5, using (1), where T max is maximal temperature in the slot, T 0 is boundary temperature on the side of the slot and Q is heat flux through the side of the slot.This detailed model can be used again for the estimation of the temperatures of conductors after the whole model is solved and the boundary conditions are known.However, this procedure is not possible for random-wound windings where the positions of the conductors are not clearly defined [5].

Finite Element Method
In order to solve this problem, the following Matlab code based on the 2D FEM was developed.To ensure that the calculation is as fast as possible, it is convenient to vectorize most of the operations.Vectorized assembly of stiffness and mass matrices is done using the code of the authors [29].To generate the triangular mesh, the third-party code from [30] is used.The result of the calculation of the temperature of the slot using this code is shown in Figure 6.

Slot-Filling Algorithm
An important part of the code is the function that fills the slot with conductors.The inputs are the dimensions of the slot, the radius of the conductors and their number N, or the required fill factor f ill.The function then returns a vector of the conductor centers.The filling loop can end in three ways: (1) the required number of conductors N R is reached; (2) the required fill factor f ill R is reached; (3) there is no space left for another conductor.
The first step is that the algorithm divides the usable area of the slot by a rectangular mesh, which can theoretically be infinitely dense, but the more points it has, the slower the calculations are.Its density is thus derived from the radius of the conductor so that the filling is relatively fast and at the same time the main aspect of the code, the probabilitydriven random selection of the centers, is not too limited.The filling loop is described in detail in Figure 7.
The selection of the conductor centers is controlled by four main variables: two logical vectors V N and V T , which contain points of the mesh belonging to two respective regions and the weights w 1 and w 2 .Vector V N contains points that are excluded from the selection of new centers, while vector V T contains points from which new centers are selected.The vector V T initially contains only a strip of points at the bottom of the slot.As more conductors are added, these vectors are updated.
The new centers are not selected from the V T vector randomly, but according to the weights w 1 and w 2 , these weights correspond to the probabilities of selection of a given point.The weight w 1 is a projection of the function f 1 (y) on the points of the rectangular mesh.This function has a maximum at the point of V T closest to the bottom of the slot.The weight w 2 is then similarly a projection of the function f 2 (x, y) on the mesh, but this time the function f 2 is a paraboloid with a maximum at the nearest point of the vector V T to the mass center of the conductors.The functions f 1 and f 2 are shown in Figure 8, this geometry corresponds to the L-type slot (Figure 9c).The orange-bounded region is V T , the red dot is then the closest point of the V T to the center of mass, which is located between the three conductors.Whether the new center will be selected according to the weight w 1 or w 2 is then controlled by the user specified parameter S d ∈ 0, 1 , which controls the density of the conductors in the slot.For S d = 0 the new center is always selected according to w 1 , for S d = 0.5 the chance is 50% for each and for S d = 1 the new center is always selected according to w 2 and the resulting distribution then corresponds to hexagonal packing with slight adjustable deviations.Two extreme cases are shown in Figure 9a,b.Filling of this slot takes approximately 0.35 s, while the size of the vectors V N and V T is 27,861 points, if we choose a denser mesh, then for 112,169 points the slot is filled in 0.95 s.

Homogenization
If we are mainly interested in the maximum temperature of the slot rather than its distribution, it may seem that homogenization is a useful tool.Its various forms, advantages and limits have already been explored by many authors, for example: [31,32].If we compare the results from the FE simulation with conductors and the simulation with homogenization for different fill factors, while keeping the same Ampere-turn value, we get the results depicted in Figures 10 and 11.These figures show the maximum slot temperatures from which the lowest maximum temperature of the whole set was subtracted for better comparison.Each of the two figures shows 5000 temperatures calculated for randomly generated conductor distributions.
The formula used to calculate homogenization is Equation (2).Where f ill is the fill factor (ratio between the area of conductors and the slot area), k a is thermal conductivity of the surrounding material and k con is thermal conductivity of the copper conductors.
These temperatures were calculated for slots designed with a fill factor of approximately 0.5, so as the fill decreases below this value, the current density increases to unrealistic values.However, this does not change the fact that the variance of the maximum temperatures increases significantly with a smaller fill factor.In the case of the L-type slot, the results of the simulation with the conductors and simulation with the homogenization start to correspond approximately from a fill factor of 0.35.The homogenization then approximates the best case, i.e., the conductors are arranged so that the maximum temperature is the lowest of all the possible arrangements.The behavior of homogenization for the slot with concentrated winding is different, here the homogenization from a fill factor of 0.4 approximate the worst case of the arrangement of the conductors.This difference in behavior is due to the fact that while in the L-type slot all the conductors are close to each other, in the slot with the concentrated winding the conductors are wound around the stator tooth and as a result there are two clusters of conductors separated by a gap.
Although homogenization can save computational time, it can be used only in cases where we are sure it matches a more accurate model with the conductors.However, this region of trust depends on several factors and will be the subject of further research.Homogenization is nevertheless used in the code to homogenize the conductor material with its thin insulating layer, which could be otherwise considered in the problem only with the use of an extremely dense mesh.

Thermal Resistances
An important function of this FE-based model is to determine the thermal resistances for the LPTN model.These resistances are calculated using (1) and the flux is then calculated from the temperature distribution according to Fourier's law (3), where k is thermal conductivity, ∇T is temperature gradient and S is a closed surface.Subsequently, the slot can be artificially divided into four parts: upper, lower, left and right, each of which passes through a part of the heat flow Q.This sub-model definition is then used in LPTN model instead of the common sub-model mentioned in Section 2.1.These thermal resistances are shown in Figure 12.
Slot thermal resistances for four directions.
The following  show the dependencies of the individual thermal resistances on the fill factor.For all resistances, a similar trend can also be seen for the temperature distribution (Figures 10 and 11), i.e., that with lower fill factors the variance of resistances increases.The strongest dependence of mean value and fill factor can be seen in case of the resistance R u -the resistance between the hot-spot and the edge of the slot adjacent to the air gap.
For resistance R d , connecting the hot-spot and the edge of the slot adjacent to the stator yoke, the dependence is relatively low.An almost linear decrease of mean value of thermal resistance with fill factor can be observed for the resistances R l and R r , which connect the hot-spot with the edges of the slot at the left and right of the stator teeth.All these values have been calculated for a slot with a concentrated winding shown in Figure 6 and the variance of values and mean curves will be different for each type of winding and slot geometry.

Experimental Validation and Measurement
To test the hybrid thermal model, we chose a traction high-speed synchronous motor with surface mounted permanent magnets on the rotor and with a randomly wound stator winding.We entered the parameters of this machine into the model and performed a thermal analysis and then measured the machine temperatures using thermo-couples.We also calculate the equivalent thermal conductivities from the machine measurements and then compare them with the calculation of the equivalent thermal conductivities calculated using conventional homogenization methods.The tested PMSM is used together with the gearbox as a compact traction machine with speed range of 3-10 kRPM.The main parameters of the motor are presented in Table 1.To make the thermal measurements as complex as possible, the machine was tested as a motor, a generator and as a stand-alone machine with its rotor removed.During the tests, the stator winding was supplied with direct current corresponding to RMS value of rated current (see Table 1) and one phase was measured.Motor was tested with a load according to duty cycle S1-continuously loaded.The thermo-couples, PT100 in four wire connections, were placed in and around the loaded phase, in several stator positions and in the stator yoke reduce noise including the shift noise.For measuring ALMEMO 5690-2M combined with ALMEMO 8690-9A were utilized.In order to accurately determine the temperature distribution in the machine, altogether 40 temperatures were continuously measured.Figure 16 shows the model of machine with some of the thermo-couples' positions marked and the Figure 17 then shows the real machine during the measurements.From these measurements it is then possible to calculate the thermal conductivities of the contact surfaces between the winding and the stator of the machine.These conductivities are calculated using (4), where l is the middle distance from the center of the slot to the contact surface, A is the area of the contact surface (area of the yoke and the area of both stator teeth) and R th is thermal resistance.These thermal conductivities are calculated for the two measured (highest 61 • C and lowest 54 • C) temperature differences between the winding in the slot and the stator tooth.Temperature differences were measured in a steady state.
These temperature differences are then considered in the calculation of the thermal resistances R th according to Equation (5), where ∆T is the temperature difference and ∆P Cu are Joule losses for rated current, in our case 625 W. The tested machine had a star-connected winding without an achievable zero node and the winding was connected in series-parallel during the tests.
Figure 18 shows the transient of 40 measured temperatures, temperatures of end winding and slot winding of loaded phase are highlighted, other values are plotted as a range of temperatures.First, the machine was powered by DC current and heated by Joule losses without a rotor.The second test was performed with the rotor inside.There was no significant difference between the measured temperatures in both cases, and it can be said that the thermal properties of the air gap are the same with or without the rotor.For the two mentioned thermal differences, the thermal conductivities are: 0.448 W•m −1 •K −1 for 61 • C difference and 0.505 W•m −1 •K −1 for 54 • C difference.

Comparison of Measured Results and Homogenization Methods
The measured thermal conductivities will now be compared with several examples of methods for calculation equivalent thermal conductivities used in the literature.Each result is calculated with the same dimensions and parameters of the machine.The first Equation ( 6) is based on the comparison of areas of individual materials in the slot.A slot is the total area of the slot, A con is area of the conductors, A ins is the area of the insulation and A imp is the area of the impregnation.Then k con , k ins and k imp are the thermal conductivities of these materials.The equivalent thermal conductivity is in this case 0.365 W The next Equation ( 7) was already presented is Section 3.2 and it is one of the most commonly used equations for homogenization.The individual parameters then have the same meaning as in the previous equation.If we want a more accurate calculation, it is possible to use the same calculation once more and calculate the equivalent thermal conductivity k con using the thermal conductivities of copper and its thin insulation.The result of this equation is 0 The next Equation (8) from [24] considers the impregnation as another insulating layer of the winding.r con is conductor radius and t ins is the thickness of the conductor insulation.The equivalent thermal conductivity is in this case 0.693 W There are many other equations and variations on these mentioned equations in the literature, but the results are then very similar.Let us now compare the measured results and the results of the LPTN model with different sub-modules for calculating the thermal resistances between the winding and the stator.
All results are compared in Table 2.The average temperature measured for a given component is always used to calculate the percentage error.The column labeled k eqv is the result of the LPTN model, into which the average equivalent thermal conductivity eqv column is then the result of the LPTN model using the equivalent conductivity of 0.512 W•m −1 •K −1 , calculated using (7), which has the smallest error of all the formulas mentioned.However, both of these calculations have an error over 10% for the temperatures of the slot winding and end winding.The calculation using thermal resistances generated for a given slot dimensions and the fill factor using the FE model has the smallest error (below 5%).In addition, after calculating the LPTN model, it allows to set the calculated temperatures as the boundary conditions and thus determine more precisely the temperature distribution in the slot, the position of the hot-spot and its approximate temperature.

Conclusions
This paper presents a new approach for determining the thermal resistances or equivalent thermal conductivities of stator slot winding.This novel method is based on a unique combination of most commonly used methods these days: LPTN and FE/CFD.The developed modular LPTN model enables fast temperature calculations in an electrical machine using a network of thermal resistors.The modularity of this model then makes it easy to redefine and add individual parts for new machine topologies.This is an important tool in the thermal design of an electrical machine, the designer can always quickly check how design changes will affect the temperatures in the machine.The FE model is used as a sub-model for accurate calculations of the thermal resistances of the slot with winding.This model is able to generate random distributions of conductors in a slot, which simulates the actual process of manufacturing of randomly wound winding.Furthermore, this model can be used to more accurately determine the temperature distribution in the slot after the boundary conditions are known from the LPTN model solution.The newly designed and implemented method allows fast thermal analyses with a high accuracy of results.The combination of these methods then achieves a margin of error below 5% compared to measurements of the temperature on a real machine.The findings from the series of measurements will be further used to further improve the software.

Figure 1 .
Figure 1.Stator slots in the cut synchronous machine.

Figure 2 .
Figure 2. Flowchart of the thermal analysis.

Figure 6 .
Figure 6.Temperature distribution in the slot.

Figure 10 .
Figure 10.Maximum temperatures of L-type slot.

Figure 11 .
Figure 11.Maximum temperatures of slot for concentric winding.

Figure 16 .
Figure 16.Position of the thermo-couples in the model.

Table 1 .
Main parameters of the measured machine.

Table 2 .
Comparison of results.