A 3D dynamic lumped parameter thermal network of air-cooled yasa axial flux permanent magnet synchronous machine

: To ﬁnd the temperature rise for high power density yokeless and segmented armature (YASA) axial ﬂux permanent magnet synchronous (AFPMSM) machines quickly and accurately, a 3D lumped parameter thermal model is developed and validated experimentally and by ﬁnite element (FE) simulations on a 4 kW YASA machine. Additionally, to get insight in the thermal transient response of the machine, the model accounts for the thermal capacitance of different machine components. The model considers the stator, bearing, and windage losses, as well as eddy current losses in the magnets on the rotors. The new contribution of this work is that the thermal model takes cooling via air channels between the magnets on the rotor discs into account. The model is parametrized with respect to the permanent magnet (PM) angle ratio, the PM thickness ratio, the air gap length, and the rotor speed. The effect of the channels is incorporated via convection equations based on many computational ﬂuid dynamics (CFD) computations. The model accuracy is validated at different values of parameters by FE simulations in both transient and steady state. The model takes less than 1 s to solve for the temperature distribution.


Introduction
YASA (yokeless and segmented armature) axial flux permanent magnet synchronous machines (AFPMSMs) have been providing challenging opportunities for industries for several years.Thanks to the absence of the yoke in the YASA type (shown in Figure 1a), it provides the highest power density and efficiency compared to other AFPMSM types, such as the toroidally wound internal stator (TORUS)-type machine shown in Figure 1b and the axial flux internal rotor (AFIR)-type machine shown in Figure 1c [1].
The proper design of this machine requires the development of accurate and fast models of the YASA machine.These include electromagnetic, mechanical, and thermal models for the machine.Therefore, accurate and fast models have been developed.Electromagnetic 3D finite element (FE) models are the most accurate but the most time consuming.Two-dimensional (2D) multi-layer FE models in [2,3] provide an alternative to the 3D FE models.However, they are still very time-consuming.Therefore, fast and accurate analytical models were developed in [4][5][6][7].These models are capable of calculating all the electromagnetic quantities (cogging torque, torque ripple, iron losses, permanent magnet (PM) losses, inductances, and no load and full load voltages) in a very efficient way.Authors in [8,9] developed mechanical analytical models for the YASA machine for a proper structural analysis of the rotor and stator for large-scale applications.In addition, the authors in [10] developed mechanical FE models for a coreless AFPMSM in a high-speed application.(axial flux internal rotor) [11].
Thermal modeling of the YASA machine is an important aspect in the design procedure.One of the most important incentives to study the thermal behaviour of electrical machines (including the YASA machine) is that the maximum torque capability of the machine depends on the maximum current density that can be reached before winding insulation failure or permanent magnets demagnetization [12].There are two types of modelling to obtain the temperature distribution inside the machine.The first type is a coupled computational fluid dynamics (CFD) and thermal FE simulation, such as in [13][14][15][16].The second type does not use CFD, and is based on appropriate equations for the convective heat coefficients, as proposed in [17][18][19][20].These coefficients can be the source inputs to a more detailed 3D FE thermal model [21] or a lumped parameter thermal network (LPTN) [22,23].
Coupled thermal and CFD analysis is highly needed with complex cooling systems.The authors in [13] studied the thermal effect of a forced air cooling using two different fan types mounted on the shaft in addition to the cooling fins mounted on the stator surface.The authors compared the cooling without fans and with fans.They proved the great capability of the thermal model to describe the temperature behaviour of the system.Additionally, Ref. [14,15] presented coupled FE and CFD analysis for studying cooling systems in electrical machines.The method was proven to be accurate.However, for the methods with coupled FE and CFD models, the computational effort is very high.
The authors in [17,18] studied the heat flow for two solid discs, in which one disc is rotating by a certain speed but the geometry is too simple for a real YASA machine.To study the air flow and heat flow in the airgaps of the YASA machine in more detail, the authors in [19] included the PMs on the rotor surface and developed analytical equations for the heat flow on the boundary surfaces taking into account the air flow inside the channels between the PMs.However, the model is not yet a complete motor model, because the solid parts of rotor and stator are not modelled.
In [21], the authors used convection coefficients of a simple geometry (rotating solid discs) using 3D finite element method (FEM) to predict the temperature in a YASA motor.However, this FEM-based model is quite slow and does not tackle complex fan-shaped magnets.
The literature lacks a full LPTN that describes the thermal behaviour of a YASA motor with innovative cooling via air channels between magnets.In this paper, the authors developed such an LPTN based on the convective heat coefficients developed in [18,19].The model is validated by FEM simulations and by experiments.Additionally, the effect of air gap length, rotor speed, and geometrical parameters of the magnet on the machine temperature is investigated by the FEM and LPTN model.

AFPMSM Lumped Parameter Thermal Modeling
This paper develops a fast parametrized 3D thermal model by coupling a separate LPTN of the stator and rotor through the introduction of the analytical convective heat transfer coefficients in [18,19].The model considers a rotor shape with holes.These holes form air channels, providing an innovative and effective way of cooling.All power losses occurring in the YASA machine are explained in Section 2.1 and are taken into account in the LPTN model.Figure 2 is a cross-section of the whole YASA machine, along with the studied rotor geometry.The stator consists of concentrated windings wound around laminated cores.The teeth are distributed around the stator circumference and housed inside a laminated aluminum housing.The two rotors are composed of a steel disc with permanent magnets mounted on the surface.Due to both the symmetry of stator and rotor construction and the thermal periodicity, only one quarter of one stator tooth and one rotor magnet pole part need to be modeled.Figure 3 shows the modeled parts of the stator and rotor.

YASA Machine Power Losses
Different electrical, magnetic, and mechanical power loss mechanisms take place in the YASA axial flux machine.These losses occur in the winding, stator core, permanent magnets, and bearings, causing a temperature increase in different machine parts.These losses are modeled by heat sources in the LPTN model.The accuracy of the predicted temperature depends on the accuracy of the power losses calculation.

Copper Losses
Copper losses Q wdg are calculated using the formulas in [4].The modeled winding is divided into three regions: the effective copper winding with power loss P wdg , upper end winding with power loss P endup , and lower end winding with power loss P enddown .The copper losses are distributed between these regions according to their respective volumes, assuming that the copper losses are equally distributed over the volume.Equation ( 1) governs the sum of the three losses parts: where N s is the number of teeth.

Core Losses
Core losses Q iron are calculated using formulas in [4].They are divided into two parts (inside the tooth P st and the tooth tips P st tips ) their sum is expressed by Equation ( 2).The core losses are also assumed to be uniformly distributed over the volume. (2)

Permanent Magnet Losses
Some losses take place in the PMs due to induced eddy currents.These losses Q PM are calculated in [6].The PM losses in one segment of the rotor equal where N m is the number of PMs in one rotor disk.

Mechanical Losses
The aerodynamic forces acting on rotary parts of the machine cause mechanical losses in the surfaces exposed to it.Extensive CFD simulations and curve fittings done in [24] have found accurate and parametrized formulas for the windage losses on every surface.These formulas are used to calculate the mechanical losses.
The equations take the form: where P f ,v and P f ,p are the windage losses for surface f due to viscous forces and pressure forces, respectively.P * f ,v , P * f ,p are the windage losses at the reference point.The following parameters are important because they influence air flow in the machine and thus the thermal behaviour: G is the gap size ratio (Equation ( 5)), Re is the rotational Reynolds number (Equation ( 5)), α m is the magnet angle ratio (Equation ( 6)), and L is the magnet thickness ratio (Equation ( 6)).Z, X are the fitting functions developed for every surface.
where Ω is the rotational speed in rad/s, D is the outer diameter, ν is the kinematic viscosity of air, s is the air gap thickness, α is the magnet angle, t is the magnet thickness shown in Figure 4, and R is the rotor outer radius.From Equation (4), the mechanical losses depend on the precise shape of the magnets and cooling channels.

Modeling of Convection Heat Transfer
The heat flux qi on each boundary surface i exposed to convective heat transfer according to Figure 3 for both stator and rotor is defined by Equation ( 7) [21].These surfaces are: stator housing outer surface, stator tooth facing rotor, rotor magnet facing stator, magnet upper part, rotor back iron facing stator, rotor back iron upper part, and rotor disc back.
where hi is the average convection heat transfer coefficient at surface i, Tsur f i is the surface temperature, and Tre f i is the reference temperature (average temperature of a nearby fluid contained within an adjacent volume i).
The reference temperature and the convective heat transfer coefficient are both found in [18,19].The reference temperature depends on the rotor and stator temperatures, and also on the parameters a i and b i : where T r is the rotor temperature, T s is the stator temperature, and T a is the ambient temperature.The convection coefficients hi besides a i and b i are dependent on the gap size ratio (G) and the rotational Reynolds number (Re) [18,19].
The convection coefficient is also a function of the Nusselt number for each surface i: where k is the thermal conductivity of air, and l i is the characteristic length.The Nu i , a i , and b i are functions of the gap size ratio and the Reynolds number.In addition to these parameters, they are also functions of the permanent magnet (PM) span ratio α m and the PM thickness ratio L for this studied rotor type.
The reference temperature is introduced in the model as a temperature source of the air near every convected surface.The convection heat transfer coefficient is used in the calculation of the thermal resistance between the air and the surface.
For surfaces with convection, the equivalent thermal resistance R c at this surface is where A c is the area exposed to the convection and h i is the convection coefficient for surface i computed by Equation (9).

Lumped Parameter Thermal Network
The LPTN is a network of thermal resistors, thermal capacitors, and heat sources connected together according to the heat flow directions and the positions of thermal nodes inside the machine.Both stator and rotor modeled segments are divided into thermal nodes, from which the heat is conducted in the axial, radial, and circumferential directions.Figure 5a is a differential representation of each modeled node in both machine segments.When there is a heat generation in the node, it is represented by the model shown in Figure 5b. Figure 5c presents the equivalent model for an element without heat generation.The central node connected to the capacitance (C) gives the mean temperature of the element.The mean temperature of the three-dimensional axis (x, y, z) gives a lower temperature than the central node.Therefore, a negative resistance equal to −1/6 of the total resistance in each direction is presented to account for this effect.This resistance comes from the independent solution of the heat conduction equation, as explained in [25,26].
The differential resistance in each direction (x, y, z) can be calculated by Equation (12).Detailed methods for the calculation of thermal resistances can be found in [27].
where dR x , dR y , dR z are the thermal resistances in the x, y, z directions.dx l , dy l , dz l are the lengths of the element.dA x , dA y , dA z are the areas perpendicular to the heat flow.K x , K y , K z are the thermal conductivities of the material.The stator core is made of thin laminated silicon steel, while the stator housing is made of aluminum laminations.These laminations are stacked together by an epoxy resin.Additionally, the stator windings are epoxy infiltrated.The thermal conductivity of the epoxy resin is low.This directly affects the overall effective thermal conductivity of the material mix, and hence the heat conduction in the machine.This material mix is taken into account by calculating an equivalent conductivity and considering the thermal conductivity with respect to the heat flow direction.
For the copper winding, the equivalent thermal conductivity K w in the lapping direction and perpendicular to the lapping direction is calculated as: in perpendicular direction, (13) where K he is the thermal conductivity of epoxy resin, K hco is the thermal conductivity of the copper, and f wi is the winding fill factor.
The thermal conductivity for both stator core and housing material mix is also calculated using Equation (13) considering the different stacking factors and material properties [21].
The thermal resistance for each geometrical node shape is calculated by integrating Equation (12) according to the direction of heat flow and dimensions.
The thermal capacitance C equals where ρ is the density of the material used, C p is the specific heat capacitance of the material, and V is the volume of the element.
In case of a material mix, the specific heat capacitance C p and the mass density ρ are calculated using where C p m is the specific heat capacity of the material mix, f wors is the winding or stacking factor, C p 1 is the specific heat capacitance of material 1, C p 2 is the specific heat capacitance of material 2, ρ m is the mass density of the material mix, ρ 1 is the mass density of material 1, and ρ 2 is the mass density of material 2. The stator segment model shown in Figure 3a is divided into two regions in the axial direction.They are at the tooth tips and at the center of the machine.Both regions are shown in Figure 6.The stator is divided into 20 subregions.The model of each subregion is constructed using the method described in Section 2.3 and Figure 5.The subregions with numbers are shown in Figure 6 with the corresponding LPTN model in Figure 7.
Figure 8 shows the LPTN for the rotor.

Solution of the LPTN
The model is solved in a time stepping manner.The temperature at each sample time is calculated using: where T n is the temperature at sample number n, t is the sampling time interval, TR ref is the reference temperature at each node divided by the resistance R d,m connected to this node, G is the thermal conductance matrix, P loss is the power losses vector, C is the thermal capacitance matrix of the machine, and T o is the temperature at sample number n − 1.
The thermal conductance matrix G is constructed as: where m is the number of nodes in each matrix, and R 1,n is the connecting resistance between each node n and node number one.
The thermal capacitance matrix equals: The power losses vector P loss : where P m corresponds to the power losses at node m if it exists.
The stator and rotor matrices are solved at each time step.However, the reference temperatures T re f used in Equation (20) are dependent on the rotor and stator temperatures, therefore, at each time step the solution is iterated with updated values of the rotor and stator temperatures until the stabilization of both.

3D FEM Thermal Model
A 3D FEM model is developed and solved.The model considers the same boundary conditions as described in Section 2. These boundary conditions are convective heat transfer for surfaces exposed to convection and thermal insulation for surfaces parallel to the cutting plane due to thermal periodicity.The heat flow equation solved by the FEM is: where ρ j is the density of material j, C p,j is the heat capacity of material j, and k j is the thermal conductivity of material j.In this FEM model, a heat generation of 1.667 W, calculated using Equation (1), is assigned to the copper part.This value corresponds to a 100 W total copper loss, which is the same value as in the experimental work.Additionally, 0.25 W of heat generation is assigned to the PM to account for a 4 W total losses.The properties of different materials used in the YASA AFPMSM machine prototype are in Table 1.The spatial distribution of the temperature of both stator and rotor parts is shown in Figure 9.This FEM model takes around 30 min to solve.

Experimental Setup and Results
Measurements were carried out using the setup shown in Figure 10 to validate the LPTN model of the YASA machine.In the setup, the 4 kW YASA machine was coupled to an induction machine of 7.5 kW rated power, 3000 rpm rated speed.The induction machine was used as a prime mover, and it was powered by a commercial drive.The induction motor speed set-points were provided to the drive by a dSPACE 1104 platform.The YASA machine geometrical parameters are presented in Table 2.The experimental setup model input parameters are in Table 3.A well-known heat source is required to validate the model, as well as its precise position; therefore, the actual Nd-Fe-B magnets were replaced by aluminum trapezoidal dummy PMs to avoid magnet losses.The windings were excited by a fixed DC current to get controlled losses and to eliminate the stator core losses.The mechanical losses are calculated using equations in [24].The mechanical losses are lumped and injected in the PM node in the rotor LPTN model.The mechanical losses are bearing losses and windage losses.They are also validated using torque measurement and dummy magnets with different shapes.Evidently, the windage losses depend on the magnet shape parameters Figure 4 (α m , L).
The stator core T core and winding T wind temperatures were measured by embedding a platinum resistance thermometer (PT100) inside them while the PM temperature T pm was measured by an infrared temperature sensor ZTP-135SR.The sensors output was sampled every 5 s for 110 min.
During the experiment, the prime mover was rotated at 1000 rpm and a fixed DC current was injected to get 100 W losses in the copper windings.The winding, core, and PM temperatures variation with time are shown in Figure 11.It can be seen that the results are in a good agreement with the experimental outcome in both transient and steady state.Table summarizes the temperatures for the winding, stator core, and PM for experimental, FEM, and LPTN model.The relative error in the winding temperature was 0.9%.On the other hand, the relative error in the core temperature was 0% and in the PM temperature it was −2.24%.It is also worth mentioning that the model took around 1 s to be solved, which is a small time that allows for rapid calculation of the machine temperatures and use in further thermal studies on the YASA AFPMSM machine.

Model Validation at Different Parameters
The model was validated at different parameters.These parameters included the air gap size ratio G, PM span ratio α m , PM thickness ratio L, and rotor speed.These comparisons give more details on the effect of these parameters on the temperature rise inside the machine.In each parameter study, all other parameters were kept fixed at the values in Table 3.

Effect of α m
Two values of α m of 0.7 and 0.9 were studied.The Results are shown in Figure 12.The LPTN is in a good agreement with the FEM results at both values of α m .The maximum relative difference between FEM and LPTN in any of the temperatures did not exceed 1.3%.The temperature has a direct relation with the PM span ratio.As the PM span ratio increases the air channel area decreases, and less air flows inside these channels, leading to lesser heat evacuation.Consequently, the temperature rises.

Effect of L
Different values of PM thickness ratio of 0.03 and 0.07 were investigated.Figure 13 shows the results at L = 0.03 and 0.07.The results are in a good agreement, with less than 1% difference in any of the temperatures.The results indicate that when L increases the temperature decreases, and vice versa.

G
The influence of the gap size ratio was investigated at 0.0067 and 0.0337.These values correspond to 0.5 mm and 2.5 mm air gap thickness, respectively.The model is in good agreement with 1.25% maximum relative difference.Figure 14 shows the results.It can be seen that the temperature increased at both values of air gap length, because the heat transfer is at some air gap length, and the heat transfer rate decreases below and above that value [18,19].

Effect of Rotor Speed
The results at rotor speed of 2000 rpm are shown in Figure 15.The results are in a good match also with 0.5% maximum relative difference.The temperature decreases when the speed increases due to higher Nusselt numbers and hence higher convection coefficients.

Conclusions
In this paper, a 3D transient LPTN model is developed for the YASA axial flux permanent magnet synchronous machine.The model considers the innovative cooling through air channels between magnets.The model takes only 1 s to solve for the temperature value in different parts of the machine: the copper winding, the stator core, PM, and other parts of the machine.The model is validated by measurements on an experimental work and by FEM results at different values of air gap thickness, rotor speed, and geometrical parameters of the PM.The results correspond very well in both transient and steady state.The model is very fast and can be used for further thermal design of the machine.

Figure 5 .
Figure 5. Node network representation.(a) Block orientation; (b) Resistance model with heat source; (c) Resistance model without heat source.

Figure 6 .Figure 7 .
Figure 6.Regions of stator modeled part.(a) Section I at the center of the tooth; (b) Section at the center of the tooth tips.

Table 3 .
Experimental model input parameters.

Table 4 .
Steady state temperatures at 1000 rpm.