*Article* **Steady State and 2D Thermal Equivalence Circuit for Winding Heads—A New Modelling Approach**

**Julien Petitgirard 1,2, Tony Piguet 1, Philippe Baucour 1,\*, Didier Chamagne 1, Eric Fouillien 2,\* and Jean-Christophe Delmare <sup>2</sup>**


Received: 22 September 2020; Accepted: 14 October 2020; Published: 18 October 2020

**Abstract:** The study concerns the winding head thermal design of electrical machines in difficult thermal environments. The new approach is adapted for all basic shapes and solves the thermal behaviour of a random wire layout. The model uses the nodal method but does not use the common homogenization method for the winding slot. The layout impact can be precisely studied to find different hotspots. To achieve this a Delaunay triangulation provides the thermal links between adjoining wires in the slot. Voronoï tessellation gives a cutting to estimate thermal conductance between adjoining wires. This thermal behaviour is simulated in cell cutting and it is simplified with the thermal bridge notion to obtain a simple solving of these thermal conductances. The boundaries are imposed on the slot borders with Dirichlet condition. Then solving with many Dirichlet conditions is described. Some results show different possible applications with rectangular and round shapes, one ore many boundaries, different limit condition values and different layouts. The model can be integrated into a larger model that represents the stator to have best results.

**Keywords:** thermal equivalence circuit; Voronoï tessellation; winding heads; nodal method; thermal resistances

#### **1. Introduction**

The study of increasingly compacted electrical machines in severe thermal environments is today an important tendency in electrical engineering [1,2]. The electrical machines with concentrated windings exhibit many advantages like high slot-filling factor, short end-winding, high fault tolerance capability, and automated winding process. Those advantages allow the high power density applications like electrical vehicles, electric aircraft, and wind turbines [3]. For such applications, accurate thermal models are necessary to describe the system behaviour. One of the main problems in the thermal study of electrical machines concerns their winding, where the temperature rises to its maximum value [2,4]. Moreover the study of thermal field becomes more and more important because the electrical designs are more compacted with more electrical density. The use of numerical tools, like the finite element method (FEM), to estimate thermal field and find hot spots in coils leads to excessive simulation time. Thus, some methods such as the nodal method like lumped thermal model have been introduced to solve quickly a thermal field in end-windings [3–8]. The objective of the present study is then to create an adaptable winding model, that reproduces a similar thermal behaviour. This model can solve all simple slot shapes with random wire layouts. Moreover, this study does not use the homogenisation technique commonly used in winding thermal calculation methods. This technique provides a homogeneous distribution of temperature while a random layout distorts this distribution

and can create other hot spots. To do this, a specific tessellation obtained via the Voronoï diagram is used. This tessellation allows evaluating thermal conductance between each adjoining wire and between a wire and its adjoining boundary. The solving is given and different applications show the results for different shapes, different layouts and different boundary conditions. The advantages of this method are to keep the fast solving from a nodal method with the exact layout of wires in winding in all simple shapes.

First, we describe four problems to be solved, the geometric choices and simplifications as well as the different simplifications applied to the model. The graph of the nodal network and the boundary conditions are given as showed with green flow-chart Figure 1. Then, in the second step, we provide two different way to estimate the thermal conductance that will be applied in the network. The first way describes a simple equation based on the shortest distance between adjoining wires. The second way gives a numerical integration which is more adapted as showed with blue flow-chart Figure 1. In a third step, the network solving is described thanks to adapted matrices. The temperature of each wire and the heat flux between each adjoining wire are solved as showed with red flow-chart Figure 1. Finally, the last step solves the thermal nodal networks of the four examples. The tool process used is described. Moreover, a comparison with Finite Volume Method (FVM) is provided to evaluate the nodal model.

**Figure 1.** Flow-chart to define different steps of the thermal model with its different possibilities.

#### **2. Thermal Equivalent Circuit**

#### *2.1. Application to a Lot of Slot Shapes with a Random Layout of Round Wires*

To the best of our knowledge, the analytic solutions for the thermal modelling of an electrical winding slot in 2D use homogenisation techniques [4–6]. These techniques use the homogeneous conductivities [9–11] for the material slot and a homogeneously distributed thermal power. A numerical solution such as the finite element method, boundary element method, finite volume method or finite difference method should be implemented. These methods refine the results and obtain a more realistic thermal distribution and also to answer the random character of wire layout in the slots. Although these solutions are accurate, they cannot meet the industrial request in terms of speed, flexibility and design purpose. For this situation, the best solving method consists of a nodal

network method. This approach is known to be faster than the others [12]. The heat sources are injected in nodes like current source in an electrical field. The difficult parts are to cut the slot into regions (stated as cells) with a homogeneous temperature and to evaluate correctly thermal resistances [13]. The thermal resistances depend on geometrical dimensions and thermal material properties. Between two nodes *i* and *j*, in steady-state, the thermal Ohm's law is applied with Δ*Tij* as temperature gradient, *Rthij* as thermal resistance and *qij* as heat flow:

$$
\Delta T\_{ij} = \mathcal{R}th\_{ij} \times q\_{ij}.\tag{1}
$$

Many geometrical and physical parameters characterise a coil. Our study concerns a section in a slot. The wires have the same geometrical and thermal characteristics. The wires are composed of insulation materials and core materials. The thin protection paper used around wires in the slot is neglected and the materials around wires can be air or resin.

A winding is composed of wires tightened to each other. However, their arrangement is not perfect. The part of wire section (∑*<sup>I</sup> <sup>i</sup>*=<sup>1</sup> *swire*,*i*) compared to the part of the slot section (*Sslot*) is obtained with the ratio in Equation (2).

$$\tau\_{slot} = \frac{\sum\_{i=1}^{I} s\_{wire,i}}{S\_{slot}}.\tag{2}$$

The ratio between slot the surface and the wire cross-section surfaces is between 0.5 and 0.8 depending on the manufacturing process [14]. The improvement of electrical machines' slot filling factors is still studied today [15]. The wire layout is hardly controlled on the machine-made end-windings. So each coil is different for the same product and the ratio is not optimal. More compact layouts are possible for the hand-made end-windings. For specific electrical machines, a specific tool can be used to obtain a flat wire layout [16]. However, this study considers a constant layout along the wire axis. Heat transfer appears only in wire layout sections. Thus the problem is reduced to a 2D study in a cross-section

If the winding is not in resin, the air is trapped between the wires. We consider during all this study that low values of 1 − *τslot* create only small air cavities. If the buoyancy forces created by the heat flow through these air cavities cannot overcome the viscous forces [17], then the air trapped in the winding is supposed to be motionless. Convection can be neglected which is ensured by a low number of Rayleigh that is, *Ra* < 1708. So, our model considers only conductive thermal transfers.

The big advantage of this model is the possibility of solving any wire layout for any slot shape, some examples are presented in Figure 2. This flexibility makes it possible to test a large number of random cases in order to detect the worst and best cases in terms of thermal heating. This gives the designer a possibility of decision for the design choices of these electric winding.

**Figure 2.** Example of different possibilities of shapes and layouts.

To create an efficient thermal nodal network, some regions are identified as uniform in temperature. Computable thermal resistances must be found between regions. Each core, made of copper, has its thermal source and its thermal conductivity *kcore* and each core is surrounded by insulation. With *kcore kins*, we will suppose that cores have a homogeneous temperature and that they will be each one a node. The temperature at the slot border in real conditions is dependant of the type of electrical machine. It should be noted that the method has been successfully modded to take into account Neumann or Robin boundary conditions. The goal of the present study is to describe the core of the method. Therefore a Dirichlet boundary condition is taken here and each boundary temperature is labelled as *Tbnd*,*i*. It should be noticed that the temperature imposed may vary spatially.

The temperature *Ti* rises due to copper losses (Joule effect) inside each core. The heat sources *Qi* from Joule losses is calculated with Equation (3) where the electrical resistivity *ρ* is supposed to be constant. *Ii* is the nominal current of each wire in steady-state and *score*, *<sup>i</sup>* the core section of each wire.

$$Q\_i = \frac{\rho \cdot l}{s\_{core,i}} \times I\_i^2. \tag{3}$$

#### *2.2. Creation of the Internal Thermal Circuit*

The heat fluxes along the wire axis are neglected. The resolution only considers the heat flow in a 2D section. This network is created thanks to Delaunay triangulation [18].

In the Delaunay diagram, each node represents a wire. All of the nodes are positioned thanks to the wire layouts. It should be noted that the determination of the layout for an industrial case is not easy and should be made via a circle packing algorithm [19]. Delaunay triangulation is applied to find all the links between the adjoining wires (blue lines in Figure 3). In the network, these links are labelled edges and they are connected to the nodes corresponding to each wire.

The nodes are the temperature potential in the corresponding wire cores and the edges could be seen as the heat fluxes between two nodes. The resolution of this thermal circuit is analogous to the resolution of an electrical circuit. So, each edge is arbitrarily oriented (blue arrow cf. Figure 3) and weighted with thermal conductances. Joules losses are added on each node like a thermal source. This thermal network is simple and its solving is very fast.

**Figure 3.** The equivalent thermal circuit for a 2D view of a rectangular slot with a random layout (**a**) and a square layout (**b**).

These thermal conductances are dependent on the materials between each temperature gradient and the thermal characteristics of these materials. The thermal conductances solving is presented in Section 3.

Some modifications of the Delaunay diagram will allow simplifying a special case and adding appropriate boundary conditions. To do this, the dual graph of Delaunay triangulation, Voronoï tessellation, is used [20] (black lines Figure 3). The Voronoï cells give several details. First of all, each node which is out of the end-winding slot is deleted and its links to each other node are cut at the slot border. Three examples are given in Figure 3a for triangles (node numbers: 10, 13, 16), (6, 9, 12) and (17, 19, 20). The semi-infinite edges of Voronoï (dotted black line) can enable to identify all the Delaunay nodes which are located at the periphery of Delaunay graph. Each peripheral Delaunay node is then connected with a new edge (green edge) to one or several boundaries. These new edges are weighted with thermal conductances adapted to the shape of the slot. In Figure 3 the green line represents the border.

Another simplifying can be done when a local square layout of Delaunay node exits, as on Figure 3b. It is considered no thermal flux between diagonal wires. So, each crossed Delaunay edge is just deleted. For example, all deleted edges are symbolized with blue dotted lines in Figure 3b.

#### *2.3. Selection of Boundary Conditions to the Thermal Circuit*

As part of this work, the boundary conditions are added directly to the inner edge of the slot. However, this nodal network can easily be added to another thermal network which solves thermal field in an electric machine. In this study, the boundary conditions are imposed at the inner edges of the slot. The thermal phenomenon on geometry borders in a thermal circuit is translated by these boundary conditions: Dirichlet condition which represents a known temperature or a Robin condition which represents a convection phenomenon between temperature fluid and walls surface. In this study, only resolutions based on Dirichlet conditions are presented. However, Robin's condition can also be used following Saulnier's recommendations [13] and an iterative resolution will have to be implemented.

A single boundary condition as presented in Figure 3 requires that all the walls of the slot are at the same known temperature. To be the most representative of the thermal environment, it is possible to add several boundary conditions. For example Figure 4 shows a possibility to apply different boundary conditions. It is possible to cut the border in another way and add as many boundary conditions as desired. It is important to check if the two boundary conditions are separated by a node of the Voronoï graph to have distinct Delaunay edges for the different boundary conditions. If not, like on each slot angle in Figure 4, the corresponding Delaunay node must be linked to the two boundary conditions in a distinct Delaunay edge. Two edges instead of one are used. A node corresponding to the boundary separation is added to the Voronoï graph. It will be used to determine the corresponding thermal conductances in Section 3.

**Figure 4.** Several boundary conditions for a rectangular slot with a random layout (**a**) and a square layout (**b**).

#### **3. Thermal Conductance Determination**

In a thermal network, the real difficulty is to find the thermal conductance between each node to estimate the heat transfer between each adjoining wire and calculate the temperature rise created by the Joule losses. This study proposes an estimation of thermal resistances which represent the inverse of thermal conductances (*Rth* = 1/*G*). It is based on the principle of the thermal bridge and provides a direct calculation of resistances.

The Voronoï tessellation is used to find the thermal conductances corresponding to each Delaunay edge (cf. Section 2). Each node is now included in a cell and each Delaunay edge can intersect with a Voronoï's cell edge. We consider 2 cases:


Figure 5a presents the first case and Figures 5b and 6 presents the second case. Also, when a cell is connected to the boundary, the Delaunay node is directly linked to the boundary and the Voronoï cell is truncated (see Figure 7). The Voronoï edge is not necessary a segment and can respect a circular slot border as shown on Figure 7c,d. In these figures, the proportion of the insulation is increased for a better understanding. It is assumed that the heat flow between two wires is only exchanged through their shared Voronoï edge.

The Delaunay nodes give uniform core temperatures of wires. The edges represent the heat flux across insulation and around media (trapped air or resin). Each material gives a thermal resistance.

**Figure 5.** Cutting an internal cell with the 2 possible cases: direct edge (**a**) and indirect edge (**b**).

**Figure 6.** Layout where wire 0 and 2 have an indirect edge.

For an internal edge (cf. Figure 5), the sum of these resistances in series provides their thermal resistances as follows:

$$Rth\_{\rm ij} = Rth\_{\rm ins, ij} + Rth\_{\rm med, j} + Rth\_{\rm med, j} + Rth\_{\rm ins, ji}.\tag{4}$$

It is assumed that the heat flow from the core to the Voronoï edge is only radial. The wires are all identical and the cells on both side of a Voronoï edge are symmetrical. So, on the same edge, the two thermal resistances of insulation are equal (i.e., *Rthins*,*ij* = *Rthins*,*ji*). The insulation resistance is given with the cylindrical known resistances [17] as follow:

$$Rth\_{\rm ins,ij} = \frac{1}{\Theta\_l k\_{\rm ins}} \cdot \ln\left(\frac{r\_{w,i}}{r\_{c,i}}\right) \,. \tag{5}$$

To determine the medium thermal resistance in a cell, it is assumed that the Voronoï edge is at a uniform temperature. With the symmetry, it is determined thermal resistances between the insulation and the Voronoï edge (*Rthmedi <sup>j</sup>* and *Rthmedji*). Two geometrical cases are possible like direct internal edge (Figure 5a) or indirect internal edges (Figure 5b). Two methods to estimate medium thermal resistances are implemented. The first method (Method A) is based on a very simple assumption which is based on the shortest path. The second method (Method B) consider that there is an infinite sum of elementary resistance in parallel between the Voronoï edge and the insulation edge.

#### *3.1. Method A: Shortest Path*

The thermal behaviour shows that the heat flow favours the easiest path. So in a uniform domain, it is the shortest path. In the *a* case, the shortest path is the distance between the two wires. So, in this segment, it is defined two identical lengths *dmini* and *dminj*. In the *b* case, the segment between two nodes does not cut the Voronoï edge, then the lengths *dmini* and *dminj* are defined by the shortest

segment between the Delaunay nodes and the Voronoï edges. It corresponds to the shortest cell border. Then it is assumed to simplify the domain with the previous cylindrical resistances which symbolize the thermal resistance between the dotted arc and the insulation arc as follow:

$$Rth\_{mcd,ij} = Rth\_{mcd,ji} = \frac{1}{\Theta\_i k\_{mcd}} \cdot \ln\left(\frac{dmin\_i}{r\_{w,i}}\right) \,. \tag{6}$$

So with symmetrical context and Equations (4)–(6), the internal thermal resistance of edge is:

$$Rth\_{cij} = 2 \times Rth\_{ins,ij} + 2 \times Rth\_{med,ij} = \frac{2}{\Theta\_i k\_{ins}} \cdot \ln\left(\frac{r\_{w,i}}{r\_{c,i}}\right) + \frac{2}{\Theta\_i k\_{mod}} \cdot \ln\left(\frac{dmin\_i}{r\_{w,i}}\right) \tag{7}$$

For boundary edges, the thermal resistances follow the same principle with some adaptations. The cell studied is defined by the Delaunay node at the wire centre and the two Voronoï nodes added at the slot limit. Moreover, the Voronoï edge can be a segment or an arc according to the slot shape. Then we have to check which is the shortest path to best determine the *dmini* as shown in the Figure 7.

When the *dmini* is found, the thermal resistances could be evaluated as follows:

$$Rth\_{c(kd-i)} = Rth\_{ins,kd-i} + Rth\_{med,kd-i} = \frac{1}{\Theta\_{l}k\_{ins}} \cdot \ln\left(\frac{r\_{w,i}}{r\_{c,i}}\right) + \frac{1}{\Theta\_{l}k\_{med}} \cdot \ln\left(\frac{dmin\_i}{r\_{w,i}}\right) \tag{8}$$

#### *3.2. Method B: Numerical Integration*

This first method gives the smallest conceivable resistance but the real value is mandatory higher. A second method is proposed to find a more precise resistance. These resistances are not directly soluble. They need a numerical integration.

The assumption made is that the flow is radial from the insulation to Voronï line. The principle is using an infinity of parallel resistances to represent the flow as shown in Figure 8. With this supposition, the effect of minimal distance is represented. Indeed, the resistance discretized for the minimum distance has the smallest value than the other radial discretized resistances. The solving of parallel discretized resistance show a value slightly higher than the smallest value.

**Figure 8.** Determination of parallel thermal resistance in a Voronoï cell.

The integral which represents this parallel resistance is written in function of angle *θ*, the origin is imposed at the centre of the study wire and the x-axis at the bottom line of the cells, as follows:

$$\frac{1}{Rth\_{\text{mech},ij}} = \int\_0^{\Theta\_i} \frac{1}{Rth\left(\theta\right)}\tag{9}$$

with

$$Rth(\theta) = \frac{\ln\left(d(\theta)/r\_{w\bar{i}}\right)}{k\_{med}\,\mathrm{d}\theta} \tag{10}$$

The distance *d* (*θ*) represents the distance between centre of studied wire and Voronï line. For internal edges or boundary edges composed by a line, this distance is written as follow:

$$d\left(\theta\right) = \frac{p}{\sin\theta - m \times \cos\theta} \tag{11}$$

*m* and *p* represent the line coefficients like *y* = *m*.*x* + *p*. If the Voronoï nodes in the local coordinate system is at this position: *Nvor*,1 (*xvor*,1, *yvor*,1) and *Nvor*,2 (*xvor*,2, *yvor*,2) so *m* = *yvor*,2−*yvor*,1/*xvor*,2−*xvor*,1 and *p* = *yvor*,1 − *m*.*xvor*,1.

For boundaries where the slot shape is circular, the distance *d* (*θ*) is written as follow:

$$d\left(\theta\right) = \max\left(r\_{cnt}, \cos\left(\theta - \phi\right) \pm 0.5\sqrt{4.R\_{bnd} - 4.r\_{cnt}^2 \sin^2\left(\theta - \phi\right)}\right). \tag{12}$$

As previously, the origin of the coordinate system is at the studied wire centre and the x-axis at the bottom line of the cells. *Rbnd* is the shape radius, *rcnt* is the distance between the origin and the centre shape position, *ϕ* is the cylindrical angle of shape centre from the coordinate system.

In this second method, the resistance of internal edges is deduced from the numerical estimation of Equation (9) with adapt application of *d* (*θ*) and the Equation (7) as follow:

$$Rth\_{\bar{c}\bar{j}} = 2 \times Rth\_{\bar{l}ms,\bar{i}} + 2 \times Rth\_{med,\bar{j}} = \frac{2}{\Theta\_{\bar{l}}k\_{\bar{l}ms}} \cdot \ln\left(\frac{r\_{w,\bar{j}}}{r\_{c,\bar{i}}}\right) + 2 \times \left(\int\_0^{\Theta\_{\bar{l}}} \frac{1}{Rth\left(\theta\right)}\right)^{-1}.\tag{13}$$

For the resistance of boundary edges with the same deduction from Equation (9) and Equation (8) as follow:

$$Rth\_{c(hd-i)} = Rth\_{ins,i} + Rth\_{med,ij} = \frac{1}{\Theta\_l k\_{ins}} \cdot \ln\left(\frac{r\_{w,i}}{r\_{c,i}}\right) + \left(\int\_0^{\Theta\_l} \frac{1}{Rth\left(\theta\right)}\right)^{-1} . \tag{14}$$

#### **4. Nodal Network Solving**

The heat sources in each electrical wire are generated by Joule losses and are noted *Qi*. This heat source *Qi* is incoming on the node *ni* representing a wire on a Delaunay graph. The temperatures corresponding to Dirichlet conditions *Tbnd*,*<sup>i</sup>* will be imposed at the boundary nodes. The temperatures at each internal node are not known and will be noted *Ti*. The heat flows through edges will be noted

*qij*. These flows *qij* are algebraic terms and can be negative or positive. In this paper we -arbitrarilyimpose the orientation of the flow from node *i* to node *j* with *j* > *i*. In this section, each linear system and matrix respects the network of Figure 4a. The nodes are listed starting with the internal nodes first then the boundary nodes like [*n*0, *n*1, ··· , *n*10, ··· , *n*20, *nbnd*0, ··· *nbnd*3]. For the edges we use a double increasing system starting with the internal edges and then the edges cutting the boundaries like [*e*0,1,*e*0,2, ··· ,*e*1,2, ··· ,*e*19,20,*e*0,*bnd*0,*e*1,*bnd*0, ···*e*20,*bnd*3].

To find the temperatures, we make an energy balance on all the nodes. The linear system obtained is composed of as many equations as there are nodes.

$$Q\_i = \sum q\_{\rm cij,out} - \sum q\_{\rm cij,in}.\tag{15}$$

The Equation (15) applied for random layout in a rectangular slot (Figure 4a) gives a linear system as follows:

$$\begin{array}{c c c c c} n\_0 & Q\_0 &=& q\_{0,1} + q\_{0,2} + q\_{0,3} + q\_{0,hd0} \\ n\_1 & Q\_1 &=& -q\_{0,1} + q\_{1,2} + q\_{1,4} + q\_{1,hd0} + q\_{1,hd2} \\ \vdots & \vdots & & \\ n\_{20} & Q\_{20} &=& -q\_{17,20} - q\_{18,20} + q\_{20,hd1} + q\_{20,hd3} \\ n\_{1nd0} & -Q\_{1nd0} &=& -q\_{0,hd0} - q\_{1,hd0} - q\_{3,hd0} \\ \vdots & \vdots & & \\ n\_{ld3} & -Q\_{ldnd3} &=& -q\_{3,hd3} - q\_{6,hd3} - q\_{9,hd3} - q\_{12,hd3} - q\_{15,hd3} - q\_{18,hd3} - q\_{20,hd3} \\ \end{array} \tag{16}$$

A specificity in thermal balance is written in boundary equation. The heat sources *Qbndi* correspond to the heat flow which leaves the system. *Qbndi* is not know but the energy conservation provide this equation: ∑ *Qi* = ∑ *Qbnd*.

If the thermal conductance is defined as *Gij* = *Rth*−<sup>1</sup> *ij* , the heat flux is developed as follows:

$$q\_{i,j} = G\_{i,j} \times \left(T\_j - T\_i\right) = G\_{i,j}.\\ T\_j - G\_{i,j}.\\ T\_i. \tag{17}$$

The linear system (16) combined with Equation (17) applied to all nodes gives a linear system where the unknowns are the temperatures *Ti*. To simplify the solving, the system is written with matrix thanks to the graph theory. First, the incidence matrix ([*Inc*] in Equation (18)) that connects edges and nodes with a sign convention. *Inci*,*<sup>j</sup>* = 1 if the heat flux leaves the node *i* and respectively *Inci*,*<sup>j</sup>* = −1 if it enters into node *i*. It should be noted that the transposed incidence matrix [*Inc*] *<sup>T</sup>* gives the two nodes connected by a specific edge. And the weighted incidence matrix [*G*] (Equation (19)) gives the thermal conductance oriented and connected at each node according to edges.

$$[Inc] = \begin{bmatrix} \varepsilon\_{0,1} & \varepsilon\_{0,2} & \cdots & \varepsilon\_{18,20} & \varepsilon\_{bnd0,0} & \cdots & \varepsilon\_{bnd3,20} \\ 1 & 1 & \cdots & 0 & \\ -1 & 0 & \cdots & 0 & \\ \vdots & \vdots & \ddots & \vdots & \\ 0 & 0 & \cdots & -1 & 0 & \cdots & 1 \\ \hline \\ 0 & 0 & \cdots & 0 & -1 & \cdots & 0 \\ \hline \\ \vdots & \vdots & \ddots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & 0 & 0 & \cdots & -1 \\ \end{bmatrix} \begin{bmatrix} 1 & & & & & \\ 1 & & & & & \\ \vdots & & \ddots & \vdots & \\ 0 & \cdots & 0 & 0 & \\ \vdots & & \ddots & \vdots & \\ \vdots & & \ddots & \vdots & \\ 0 & \cdots & 0 & 0 & \\ \end{bmatrix} \tag{18}$$

$$[\mathbf{G}] = \begin{bmatrix} e\_{0,1} & e\_{0,2} & \cdots & e\_{18,20} & e\_{\mathrm{fmd},0} & \cdots & e\_{\mathrm{fmd},20} \\ \hline G\_{0,1} & G\_{0,2} & \cdots & 0 & \mid & G\_{\mathrm{fmd},0,0} & \cdots & 0 \\ -G\_{0,1} & 0 & \cdots & 0 & \mid & G\_{\mathrm{fmd},0,0} & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots & \vdots & \mid & G\_{\mathrm{fmd},0,0} & \cdots & 0 \\ 0 & 0 & \cdots & -G\_{\mathrm{fmd},18,20} & 0 & \cdots & \mathbf{G}\_{\mathrm{fmd},2,0} & \mathbf{0} \\ \hline 0 & 0 & \cdots & 0 & -\mathbf{G}\_{\mathrm{fmd},0,0} & \cdots & 0 & \vert & \mathbf{n}\_{\mathrm{fmd}} \\ \vdots & \vdots & \ddots & \vdots & \vert & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & 0 & \vert & 0 & \cdots & -\mathbf{G}\_{\mathrm{fmd},20} \end{bmatrix} \begin{matrix} n\_0 \\ n\_1 \\ \vdots \\ n\_{\mathrm{f}} \\ \hline \end{matrix} \\ \tag{19}$$

With these different matrices which come from the graph theory, the vector of edges temperature potential is given with Equation (20). This vector coupled with matrix [*G*] lets us write the previous linear system to matrix system as Equation (21) with an efficient and computable tool:

$$\mathbb{E}\left(\left[Inc\right]^T \cdot \left[T\right]\right)^T = \begin{bmatrix} \varepsilon\_{0,1} & \varepsilon\_{0,2} & \cdots & \varepsilon\_{Ind0,0} & \cdots & \varepsilon\_{Ind3,20} \\ T\_0 - T\_1 & T\_0 - T\_2 & \cdots & T\_0 - T\_{\rm{bind0}} & \cdots & T\_{20} - T\_{\rm{bind3}} \end{bmatrix} \tag{20}$$

$$\underbrace{\left[\mathcal{G}\right] \cdot \left[Inc\right]^{T}}\_{\left[M\right]} \cdot \left[T\right] = \left[\mathcal{Q}\right].\tag{21}$$

This matrix system is detailed in Equation (22). It can not be solved directly. In the vector [*Q*], the Joule losses are known but the heat flows out of the system on boundary nodes are unknown (*Qbndi*). However, the number of unknowns (temperature of the internal nodes) corresponds to the number of internal nodes. It is not mandatory to keep the corresponding equations at the boundary nodes to solve. So, the first simplification, we remove equations from boundary nodes in [*M*] and [*Q*] (i.e., equations from *nbnd*<sup>0</sup> to *nbnd*<sup>3</sup> in grey part). Then the system is horizontally split into 2 matrices labelled *ML*<sup>∗</sup> and *MR*∗. This allows to separate the unknown internal temperatures (red part) from the unknown boundary condition temperatures (green part).

Equation (22) could be rewritten as:

$$M^{L\*} \times T^L + M^{R\*} \times T^R = Q^\*. \tag{23}$$

Equation (23) could be easily transform a now solvable Equation (24) with *ML*∗ , [*Q*∗], *MR*∗ , *TR* known and *TL* unknown.

To determine all the heat fluxes in the system ! *qedge*" we could rely on the vector [*T*] previously determined. We just create the conductance vector ! *Gedge*"

$$\mathbf{a}\begin{bmatrix} \mathbf{G}\_{\text{edge}} \end{bmatrix}^T = \begin{bmatrix} \mathbf{G}\_{0,1\prime} \, \mathbf{G}\_{0,2\prime} \cdots \, \mathbf{G}\_{1,2\prime} \cdots \, \mathbf{G}\_{18-20\prime} \, \mathbf{G}\_{\text{bmd}0,0\prime} \cdots \, \mathbf{G}\_{\text{bmd}3,20} \end{bmatrix} \tag{25}$$

with the Equations (17), (20) and the known vector [*T*], we could determine ! *qedge*" thanks to the following equation:

$$
\begin{bmatrix} q\_{cd\,\text{g}\epsilon} \end{bmatrix} = \begin{bmatrix} \text{Inc} \end{bmatrix}^T \cdot \begin{bmatrix} T \end{bmatrix} \cdot \begin{bmatrix} \text{G}\_{cd\,\text{g}\epsilon} \end{bmatrix} \tag{26}
$$

#### **5. Application and Validation**

Four applications corresponding to Figure 2 with the two thermal conductance determinations (cf. Section 3) are given in this section. The thermal properties, slot dimensions and wire positions are mentioned in Appendix A. These results show different temperature fields with different coil implementations. These cases are not representative of existing electric coils in an electrical machine in terms of dimensions but they tend to prove the ability of the method to be applied to several applications.

Figure 9 shows a clear understanding of tool process presented in this paragraph. A comparison with a Finite Volume Method (FVM) used in ANSYS FLUENT is done with the same boundaries and assumptions. All materials are solid domains, so only the thermal equation is solved. The domain discretization uses triangle elements and it is applied with GMSH meshing [21]. The common interface between the different materials is meshed with conforming mesh. The mesh is refined at core and insulation boundaries until the results are independent of the mesh. For the nodal model, to apply the mathematical process, the PYTHON code is used with some libraries. First, the random layouts are generated from a 2-dimensional real-time rigid body physics engine: PYMUNK library as described in a previous work [22]. The Delaunay and Voronoï networks are created thanks to Nocaj study [20] with different library tools in SCIPY library like Convexhull. The NETWORKX library [23] is used to save and transform Delaunay and Voronoï network. Each node and each edge is provided with full data like position, weight, dual edges id ... These tools provide the incidence matrix and the weighted incidence matrix. The Quadpack routine [24] via the SCIPY library is used to solve numerical integrals in the determination of thermal resistances. Finally, Equation (24) is solved with a standard linear algebra routine from Scipy based on the GESV Lapack routine [25].

**Figure 9.** Software, programs and codes used to implement models, methods and comparison.

All the contours of the nodal model results are obtained with a radial basis function (Rbf) interpolation [26] on a sufficiently refined grid. The Rbf interpolation is directly available in SCIPY library and could be easily plotted via the MATPLOTLIB library [27].

To guide the interpolation, all the node's temperatures are known as the boundary condition ones. Also, we could determine several temperatures on the insulation of each node. As an example to determine the temperature of the insulation of the node *i* along the Delaunay edge between node *i* and *j* (*Tins*,*i*,*j*), we could use the heat flux *qij* and *Rins*,*i*.

$$T\_{\rm ins,i,j} = T\_i + R\_{\rm ins,i} \times q\_{ij}.\tag{27}$$

These additional temperatures are used in the interpolation process to obtain a more detailed contour, as shown on Figure 10 and in Appendix B.

#### *5.1. Comparison between Two Conductance Methods and a Commercial Software*

To identify the different results the indices *A* and *B* correspond to the model result with respectively the first method and the second method cited in the section of conductance determination. The acronym *FVM* corresponds to commercial software results. In all this section, each heat-up and each relative difference are based on the temperature: *Tbnd*0.

Table 1 shows that the model A with thermal conductances based on the minimum length underestimates the temperature compared to the results of the FVM. This solution is not protective and the gaps are large. The relative gaps are around 30.8% on the random layout and 18.3% on the square layout. These gaps are up to 30 ◦C on a heating estimated by the FVM of 137 ◦C (wire 10 of the square layout). Despite these gaps, the qualitative distribution of the temperature is respected (see contour Figure A1 and Figure A4).

**Table 1.** Results for square and random layout in rectangular slot with the model A, B and Finite Volume Method (FVM) compared. Corresponding thermal contour in Appendix B respectively in Figure A1 and Figure A2.


Model B always overestimates the temperature. For the square layout the deviation is approximately 13.3% and for the random layout is around 12.0%. In the random case this gaps ranging from 4.8 ◦C (wire 19) to 15.6 ◦C (wire 9) for heat-ups ranging from 55.7 ◦C (wire 19) to 121.1 ◦C (wire 11).

The observations for the round slots are similar (Table 2). With method A, the temperatures are underestimated. The relative gaps for the hexagonal layout are around 63.5% and the relative gaps for the random layout are around 42.5%. In the random case, the heat-up differences between the method A and FVM range between 26.4 ◦C (wire 17) and 61.5 ◦C (wire 6). The FVM corresponding heat-ups range between 64.7 ◦C (wire 8) and 145.9 ◦C (wire 6).

The B model is still protective for this slot shape. For the hexagonal layout, the relative deviation is around 6.5%. The heat-up differences range between 3 and 4.5 ◦C for heat-up between 43.3 and 79.7 ◦C. For the random layout, the relative gap is around 12.1% with heat-up differences of 8.0 ◦C (wire 8) to 16.4 ◦C (wire 6) for FVM heat-up between 64.7 ◦C (wire 8) and 145.9 ◦C (wire 6).

Method A has the advantage of being very fast with a direct calculation of thermal conductances. It qualitatively represents the thermal behaviour but greatly underestimates the temperatures. Method B is always protective, more precise and the temperature distribution is also preserved (Appendix B). Its disadvantage is the estimation of the integral when calculating the thermal conductance. This solution requires a bit more IT resources but remains faster than the solving of finite element methods.


**Table 2.** Results for hexagonal and random layout in round slot with the model A, B and FVM compared. Corresponding thermal contour in Appendix B respectively in Figure A3 and Figure A4.

#### *5.2. Comparison between the Different Shapes and the Wire Layouts*

For the next results, only the method B is presented and discussed. For the Figure 10a,b four distinct boundary conditions are applied such as *Tbnd*,0 = 50 ◦C, *Tbnd*,1 = 80 ◦C, *Tbnd*,2 = 65 ◦C and *Tbnd*,3 = 62 ◦C. For the Figure 10c,d, only one boundary temperature is imposed at *Tbnd*,0 = 50 ◦C. For all cases in this figure, the electrical current in each wire is 7.5 A which corresponds to a heat flow of 13.15 W to Joule losses.

**Figure 10.** Temperatures contour for 2 shapes and 2 layouts. Results data from Tables 1 and 2.

The results prove that the modelling technique can provide individual wire temperatures. With a central area hotter than the periphery, the thermal field between (a) and (b) are very similar: the hotspot is in the centre of rectangle, that is, on wires 7, 10, 13 for the square layout and on wires 5, 8, 9, 11, 13 and 14 for the random layout. However, in detail, the number of wires in the hot spot is 3 for a square layout against 6 for the random layout. So, in the random layout, the temperature is lower (180 to 185 ◦C against 202 to 207 ◦C) but the hot spot is spread over a larger area. The boundary conditions greatly influence the temperature field. Consequently, in a random geometry, a bigger distance to the border decreases the influence of the boundary temperatures. This influence is visible between wires 18, 19, 20 of Figure (a) and wires 0, 1, 3 of Figure (b). Indeed in Figure b wire 0 is very close and has a minimal temperature influenced by the boundary conditions while wires 1 and 3 are warmer. The same phenomenon is visible in Figures (c) and d when we consider wires 15, 17, 18 to (c) and 0, 1, 2, 3 to (d).

The heat-up between the hexagonal (c) and random (d) layout for the round shapes is not of the same order of magnitude. The heat-up in the hexagonal layout is very lower than the random layout. Indeed the hexagonal layout optimises a small and constant space between each wire. The insulation created by the air is minimised. Many wires are close to boundary and they help the heat flow to escape.

The comparison between the shapes is interesting despite the lower number of wires in the round geometry (19 wires) than in the rectangular geometry (21 wires). The comparison shows that the round shape has the lowest (c) and the biggest (d) global temperature between all the cases. This shows the importance of the wire layout and the choice of shapes to design the coil.

#### **6. Conclusions**

This study provides a methodology for analysing the heat-up of any set of wires in a end-winding or any other device. Based on geometrical assumptions (Delaunay triangulation and Voronoï tessellation), the model creates a thermal network that could solve easily the temperature field. The evaluation of the thermal transfer thanks to thermal resistances and conductances is also provided. A second estimation more precise of thermal resistances and conductances is provided. Network set-up, adaptation, matrix writing and resolution are detailed. The model has been compared with FVM and several experiments are planned. To refine the model several issues should be tackled:


**Author Contributions:** Conceptualization, E.F. and P.B.; methodology, J.P. and T.P.; software, J.P., T.P. and P.B.; validation, D.C.; investigation, J.P. and T.P.; data curation, T.P.; writing—original draft preparation, J.P.; writing—review and editing, P.B., D.C. and E.F.; visualization, D.C.; supervision, D.C. and E.F.; project administration, J.-C.D. and D.C.; funding acquisition, J.-C.D. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Acknowledgments:** The authors express their gratitude to the Groupe PSA, to the EIPHI University Research School (contract "ANR-17-EURE-0002") and to the National Association for Research and Technology (ANRT convention "2017/1091") for their support in carrying out this work.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A. Application: Materials Properties, Dimensions and Wire Positions**

*Appendix A.1. Materials Properties*


#### *Appendix A.2. Slot Dimension, Slot Properties and Wire Positions*

The Table A1 gives all geometry data of wires and slots corresponding at rectangular and round slot.


**Table A1.** Type and dimension of slot application.

The Table A2 gives all wire positions for the 4 slot applications.

**Table A2.** Wire positions in 4 slot applications cf. Figures 2 and 10, the origin point is at the slot gravity center.


#### **Appendix B. Results Data**

The Figures A1 and A2 give respectively a temperature comparison between the two methods and FVM for the rectangular shape with a square layout and the rectangular shape with a random layout. The Figures A3 and A4 give respectively a temperature comparison between the two methods and FVM for the round shape with a hexagonal layout and the round shape with a random layout.

**Figure A1.** Comparison of thermal contours for rectangular shape with a square layout between method A (**a**), method B (**b**) and FVM (**c**). Results data from Table 1.

**Figure A2.** Comparison of thermal contours for rectangular shape with a random layout between method A (**a**), method B (**b**) and FVM (**c**). Results data from Table 1.

**Figure A3.** Comparison of thermal contours for round shape with a hexagonal layout between method A (**a**), method B (**b**) and FVM (**c**). Results data from Table 2.

**Figure A4.** Comparison of thermal contours for round shape with a random layout between method A (**a**), method B (**b**) and FVM (**c**). Results data from Table 2.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
