Development of a Thermal Analysis Tool for Linear Machines

The thermal design is one of the most important stages in the design process of electrical machines. Thanks to software packages, like Motor-CAD, rotating machine designers can predict the thermal operation of the machines with high precision. However, a Motor-CAD equivalent for linear machines does not exist. Thus, linear machine designers must develop specific thermal analysis tools when designing the machines. In this article, a generic thermal analysis tool for different kinds of linear machines is presented. The model has been designed in MATLAB Simulink. Hence, it should be easy to implement for most engineers. The article describes the configuration of the different elements of the tool. The calibration parameters and procedure, and typical values of the calibration variables, are also given in the document. Finally, in order to demonstrate the generic nature of the tool, the model is experimentally validated via DC thermal tests to a linear induction machine and a linear switched-flux permanent magnet machine. The results show that, despite being simple and easy to implement, the model can predict the thermal operation of different machines with high precision.


Introduction
In the last few years, thermal optimization has become a really important stage in electrical machine design. The thermal sizing of electrical machines was previously based on several rules of thumb and rudimentary variables. If the maximum allowable temperatures were surpassed, the designers would increase the size of the machines to enhance heat dissipation. Hence, increased cost and volume had to be accepted as the only solution to overheating issues [1].
Thanks to the latest thermal analysis tools, the designers are able to make precise predictions of the working temperatures at different parts of the machines. Thus, both the cost and thermal operating point of the machines can also be optimized in the design stage. There are currently two main types of tools for the thermal analysis of electrical machines, numerical and analytical. The numerical tools divide the analyzed geometries into small samples. Then, the behavior of each sample is obtained with the physical formulations that describe the phenomenons that are taking place in each of them. These methods allow the analysis of complex geometries. However, they are quite demanding to set up, and the required computation time is significant, especially when complex geometries are analyzed [1][2][3][4].
Computational Fluid Dynamics (CFD) techniques are one of the most common numerical methods for the thermal analysis of electrical machines. In this case, the tools compute fluid flows and heat transfer between solids and fluids [5]. It is hard to calculate the heat transfer coefficient analytically in regions with complex geometries, such as the end windings or the airgap. Thus, CFD simulations are often used to predict the operating temperatures of electrical machines [5][6][7][8][9][10][11]. However, as this technique requires high com-the necessary guidelines for linear machine designers to implement the generic model, in order to use it as a reference when designing the machines. In order to validate the generic nature of the tool, the model is validated via several thermal tests to a linear induction machine and a linear switched-flux machine. Some recommendations are also given for the calibration procedure of the thermal model.
The article is organized as follows: Section 2 introduces an overview of the thermal model. Section 3 describes the thermal networks of the different elements of the machine. The validation of the model is given in Section 4. The results are discussed in Section 5, and, finally, the conclusions of the work are listed in Section 6.

Overview of the Thermal Model
Simulink contains a toolbox that is prepared for LP thermal network modeling. Thus, it is a very appropriate software package for a fast implementation of this kind of model. Figure 1 shows the general view of the implemented network.  The toolbox has predefined elements for the three main heat transfer mechanisms: conduction, convection and radiation, and also thermal capacitances. The conductive thermal resistance R cond is calculated with where l, k, and A are, respectively, the thickness of the material in the direction of the heat-flow, the thermal conductivity of the material and the area of the plane of the material that is perpendicular to the heat-flow. In Simulink, the thermal resistance is calculated by the thermal resistor block itself, and the user must specify the values of the thickness, the area and the thermal conductivity. However, the block does not allow to change the thermal conductivity of an element in function of the temperature of the material. Thus, certain elements have been modeled with variable thermal resistors instead of bare conductive thermal resistances. As an example of the variation of the thermal conductivity of materials, Figure 2 displays the evolution of the thermal conductivity of air in function of its temperature.  All the thermal resistances that model some kind of air volume have been consequently modeled with the aforementioned variable resistors.
The same issue has been detected with the convective heat transfer blocks of the library. They only allow a single value for the heat transfer coefficient. Thus, variable resistors have also been used for modeling all the convective heat-flow paths. The convective resistance R conv is calculated with where h is the convective heat transfer coefficient, and A is the total surface area of the body in contact with the surrounding fluid. The value of h is dependant of the shape, size, and temperature of the analyzed object and also the thermal properties of the surrounding fluid, air in this case. Thus, it is not possible to use a single value to model the thermal behavior of the machine precisely enough. Regarding radiative heat transfer, the variation of the heat-flow is only driven by the temperature difference of the elements in the heat-exchange. Thus, the pre-defined block from the library of Simulink is valid. The block establishes a heat flow Q based on Here, A, T s , and T ∞ are the emitting area, the temperature of the emitting body and the temperature of the surrounding surfaces. k rad is the radiation coefficient. This last parameter is also changed by the shape of the objects. However, the formula for parallel plates, has been used for all the radiating elements to simplify the calculation process. In this way, the only configuration parameters are the area of the body and the emissivities of the emitting ε 1 and receiving materials ε 2 . σ is the Stefan-Boltzmann constant, 5.6704 × 10 −8 W /m 2 K 4 .
To ensure universality between different devices, the model only consists of three parts: end windings, primary core, and slot. The secondary is not contemplated by the model because, as demonstrated in Reference [22], the temperature of the secondary side of the machine is not affected by the losses if the stroke length is large enough.

Thermal Network of the Slot
The equivalent thermal network of the slot from Figure 3 is based on the layered winding thermal model of Motor-CAD [27]. The advantage of this model is that it is precise enough to detect hot spot problems in the windings, but the calibration of the resistance is only dependant on the thermal conductivity and the quality of the impregnation, and the fill factor of the slot.  The thickness of the copper+air+impregnation layers is the same as the diameter of the conductors of the slots. In the initialization process, the the thermal analysis tool starts from setting a low thickness for the air+impregnation layers, and this value is iteratively increased until the desired fill factor is obtained.
Each element of the slot has its own thermal resistor and capacitor as it can be seen in Figure 3b. In order to limit the height of Figure 3b, the thermal network of a machine which has a single copper and air+impregnation layer is shown in the figure.
Thanks to the winding model, it is also easy to define the values of the conductive thermal resistances. The thickness of the layers is taken as l, and A is calculated by where h layer and w layer are the average height and average width of the layers.
In the case of the air+impregnation layers, the conductive thermal resistance depends on the quality of the impregnation q imp . This parameter contemplates the small air gaps that can appear in the impregnation due to imperfections in the manufacturing process of the machine. The value of the conductive resistance of this layer R AirImp is calculated following where k imp and k air are the thermal conductivities of the impregnation and the air.
The thermal capacitances C are calculated from the mass and the specific heat capacity c p of the materials, Here, ρ is the density of the material of the layer, and V the total volume of the layer. The specific heat capacities, thermal conductivities, and densities of some typical materials are given in Table 1. The thermal properties of air at 1 atm are given in Appendix A.

Thermal Network of the Primary Core
As it can be seen in Tables 1 and A1, the thermal conductivity of the primary core (steel) is much higher than that of the impregnation and the air. Thus, the conductive thermal resistance of the primary core has been neglected. This means that the model assumes that the temperature of the core of the machine is the same along the whole primary. Due to this simplification, the temperature of the copper is also assumed to be the same in all the slots of the same type because they are exposed to the same losses and the same thermal resistances.
The same happens with the influence of the end effect. In the reality, the amplitude of the airgap flux density is different in the edges of linear machines when compared to that in their central part. This generally creates a non-uniform distribution in the core losses. However, the influence of these non-uniform distribution of losses over the temperature is generally small. Hence, assuming a uniform temperature in the whole core of the machine, is a good approximation that helps to the simplification of the LP model.
The thermal network of the primary, which is shown in Figure 4, consists of the thermal capacity of the primary core, and convective and radiative thermal resistances to model the heat-exchange with the surrounding elements. The primary core has been separated in different elements to allow an easier calculation of the convective heat transfer coefficient and the heat exchange areas. These elements are the top and bottom surfaces, the left and right sides, and the front and rear faces. Figure 1a shows the simplified geometry that has been assumed in the construction of the model.  The first step when calculating the convective heat transfer coefficient is to determine if it is a natural convection environment or a forced convection situation. In the examples of this article, the thermal tests for the validation of the model were made statically. Hence, natural convection correlations have been applied.
Then, the thermal properties of the fluid are calculated at operating temperature, The temperature for the calculation of the properties of the fluid is T f , and T s and T ∞ are the temperatures of the surface and the ambient in Kelvin degrees. As both the ambient temperature and the surface temperature can change while the machine is operating, in the developed tool, the user has to specify the range of ambient and surface temperatures for the simulation. The tool then generates a matrix of T f for every possible combination. Then, the respective matrices are calculated for the properties of the fluid, based on Table A1.
Once that the properties of the surrounding fluid are obtained, the next step is to calculate the thermal volumetric expansion coefficient of the fluid β. For ideal gases, this follows Then, the Grashof number Gr L and the Rayleigh Ra L number are obtained from and where g is the gravity in m/s 2 , L c is the characteristic length of the surface, v is the kinematic viscosity of the fluid at the operating temperature, and Pr is the Prandtl number. Later, the Nusselt number is calculated from empirical correlations. In this case, the correlations that have been used are those presented in Reference [26] for vertical and horizontal plates. The correlations are given in Table 2.
In the final result, h is a matrix with one value for each T s and T ∞ combination. This h is introduced to the model in lookup tables, and the model interpolates the value of the convective heat transfer coefficient depending on the ambient temperature and the temperature of each surface.
The emissivity of the materials has been extracted from Table A-18 in the appendix of Reference [26]. The given emissivity values are dependant on the finishing of the material; thus, approximated values have been selected. These values are collected in Table 3.

Thermal Network of the End Windings
Similar to the thermal network of the primary core, the thermal conductivity of the end windings has been taken as infinite. Thus, the temperature of the copper of the end windings is the same as the average temperature of the copper in the slot. The only thermal resistances in the thermal network of the end windings ( Figure 5) are those corresponding to the radiation and the convection. This is a simplification that has been made to reduce the complexity of the model. Obviously, there is an error in the calculation of the temperature of the copper in the end windings. But this can be assumed to be generally small. However, when the heat flow through the primary core is very high, the end windings can contain the actual hot-spot of the copper. Due to the simplification, this kind of hot-spot issue cannot be detected with this model. As it can be seen in Figure 1a, the end windings have been modeled as rectangular prisms. The dimensions of the prism depend on the type of winding that is being modeled. Figure 6 shows the length of the prism (see l ew in Figure 1a) for a double-layer concentrated winding and a single layer distributed winding. The height of the prism is the height of the slot minus the thickness of the wedge, and the length is obtained by subtracting the width of a primary tooth to the active length of the machine. The convective heat transfer coefficients for the faces of the prisms are calculated with the process from Section 3.2.

Examples of Layouts
This tool can be used to analyze the thermal behavior of linear machines in different layouts. In this section, some examples of layouts are given, and the correlations that should be adopted in each of them are listed. An aluminium sheet has been introduced to the model in order to contemplate the mechanical system which the machines are attached to. The layouts are depicted in Figure 7.
When the layout is changed, the convective heat transfer coefficient matrix must be recalculated with the correct correlations for the defined configuration. Tables 4-6 show the correlations for the layouts that are displayed in Figure 7.

Validation of the Thermal Model
Two different linear machines have been used to validate the thermal analysis tool: a linear induction machine, and a linear switched-flux machine. The idea behind this double validation was to demonstrate the generic nature of the tool. The main parameters of the machines can be found in Table 7. Figure 8 shows the prototypes that were tested. The LIM has a double-layer distributed winding, and 3 slots at each side of the machine are half-filled. On the other hand, the LSFPM has a double-layer concentrated winding. The width of the edge slots is half of that of the rest of slots. The end slots are opened, so there are no teeth at the ends of the machine. The validation of the thermal model was done via DC thermal tests. In these tests, the machine was supplied with a constant DC power, with the phase windings connected in series. This way, the magnetic losses and mechanical friction losses were avoided, and the overall heat flow of the machine was controlled directly with a power supply. The main advantage of DC thermal tests is that the only voltage in the phase coils is generated because of the resistance of the armature winding. Thus, if the current and the voltage are measured, the value of the resistance can be calculated with Ohms law, Then, the resistance can be used to obtain the average temperature of the copper of the machine. The resistance of the phase windings varies according to where R 0 is the resistance of the series connected phase windings, α cu is the temperature coefficient of the copper, 0.00393, and T 0 is the temperature at which R 0 was measured.
Assuming that the machine is at ambient temperature at the beginning of the tests, T 0 can be considered to be the same as the ambient temperature. It can be deduced from (14) that the first step when testing the machines is to measure the initial resistance of the phase windings and the ambient temperature when the machine is cold. Then, the data from the tests can be used to calculate the temperature of the copper with Thermocouples were used in the rest of the points where the temperature was wanted to be known. The placement of the thermocouples in the thermal tests to the LIM and the LSFPM is shown in Figure 9. And the temperatures of the machines under different conditions are given in Figures 10 and 11. In the case of the LIM, the heat dissipation capability was different depending on the layout. At 135 W, the copper of the machine reached 120°C in the horizontal test. However, the same heat was dissipated at a lower temperature when the machine was mounted on an aluminium plate. In fact, the heat dissipation capability in the tests was 44.4% higher with the machine mounted in the aluminium frame. In this configuration, 195 W were dissipated with the copper at 120°C. This is consistent with what was said by Staton et al. in Reference [15,16]. Av   TC01  TC02  TC03  TC04  TC05  TC06  T amb   T cu Av   TC01  TC02  TC03  TC04  TC05  TC06  T TC01  TC02  TC03  TC04  TC05  TC06  T amb   T cu Av   TC01  TC02  TC03  TC04  TC05  TC06  T   The aluminium sheet also helped to reduce the temperature of the core of the machine. According to assuming that the total slot-to-core thermal resistance R th remains equal, the temperature difference ∆T between the coper and the core is higher when there is a higher heat flow Q through the slot. The test platform of the LSFPM can be observed in Figure 12. The coil holders of this machine were 3D printed with RS PRO's PLA filament. This material has a vicat softening temperature of 60°C [30]. Therefore, the maximal temperature of the copper was kept under 85°C to ensure the integrity of the machine, and avoid an excessive deformation of the holders. Two main temperature groups can be observed in Figures 10 and 11. In order to compare the test results to the model, the temperatures of those groups, copper and core/magnets, have been averaged.
The model of each machine was calibrated by adjusting the impregnation quality and the convective heat transfer coefficients. The results of the calibrated models are compared to the results of the tests in Figures 13 and 14.

Discussion
It is clear from Figures 13 and 14 that the thermal model that is described in this article is able to predict the thermal performance of different types of linear machines precisely. The model can also predict the thermal performance of the machines in different layouts.
Although the model calculates the values of the thermal resistors based on the geometric parameters of the machine, some parameters must be calibrated.
The main calibration parameters are the quality of the impregnation, and the convective resistances of the core and the end winding. The value of the impregnation quality has to be coherent with the fabrication process of the machine. Boglietti et al. analyze the parameter sensitivity of the thermal response of some totally-enclosed fan-cooled induction machines in Reference [4]. For the impregnation goodness factor, the authors report values between 0.4 and 0.5 in machines where the impregnation process is not optimized. Thus, this value is a good starting point when there is no previous reference about the impregnation quality of the employed process. If more advanced techniques, such as vacuum impregnation, are adopted, the impregnation quality can be enhanced up to a value of 1 [15]. When calibrating the model, the impregnation quality is used to define the difference between the average temperatures of the copper and the primary core. However, a higher goodness factor will also increase the temperature of the primary core. Thus, if the temperature of the primary core gets too high, the convective heat transfer of the end windings has to be increased to reduce the heat flow through the core.
The modification of the convective heat transfer coefficients of the core is another calibration alternative. If the temperature of the core is too high, there are two ways to lower its value. The first one is to enhance the convective heat transfer coefficient of the primary core. This will increase the thermal jump from the copper to the core. If the thermal jump gets too high, then the end winding convection must be corrected.
If the resultant correction factors for the convective heat transfer are too high or too low, the impregnation quality should be reviewed. As a reference, Boglietti et al. [4] apply correction factors in a range of 0.87 to 1.35 to the convective heat transfer of the models of five different induction motors.
There is no such thing as a generic mechanical system. The mechanisms are specific to each machine and its application. Therefore, the model cannot predict the influence of the external structure of the machine. In this case, the increased heat-flow due to the aluminium plate of the LIM and the mechanical system of the LSFPM are modeled as a reduced convective thermal resistance in the surface of the yoke.

Conclusions
In this article, a generic thermal analysis tool was developed for linear machines in MATLAB Simulink. The model is based on simple shapes and a layered winding model to allow an easy calculation of the thermal resistance and capacitance values.
The article delivered the necessary guidelines for an easy implementation of the model in Simulink. Moreover, the calibration procedure and typical values of the calibration parameters were given. Thus, this article should serve as a reference manual for linear machine designers when implementing this kind of tool.
The model was experimentally validated with a linear induction machine and a linear switched flux machine, and the calibration procedure of the model was also explained. The results show that the model is precise for both machines, and in a handful of different scenarios. Therefore, the generic model of this tool was validated.
Simulink has been demonstrated to be a very useful tool when developing this kind of tool. The predefined blocks and the automatic solving make it very easy to implement this model. However, the convective heat transfer and the conductive thermal resistance blocks do not allow adaptive parameters. Therefore, these elements must be modeled as variable thermal resistances.
The influence of the mechanical system is not contemplated by the model. Thus, when designing a new machine, the authors would recommend to build a dummy prototype with an approximated size, and characterize the mechanical system. In this way, an accurate prediction of the thermal performance is ensured from the design stage.