**By R. Belmans**

**Coupled thermo-electromagnetic problems** have to be considered in the simulation of realistic electromechanical devices and electroheat installations. To handle this type of problem, an interface program to coordinate the finite element method simulations and to perform intermediate calculations, such as heat source evaluation, numerical relaxation and mesh transitions, is recommended.

**The material data used to define** FEM problems are often strongly dependent on the temperature. Since material data parameters occur in many coefficients of the electromagnetic field equations, the calculation of the electrical and/or the magnetic field coupled to the thermal field is recommended. The heat generation of electromagnetic nature results in a coupling of the source term of the right-hand side of the thermal equation. The combined problem in the electromagnetic and the thermal domain is generally described by Helmholtz-like differential equations. The discretisation yields two or more sets of algebraic equations that have to be coupled numerically: the electric field and/or magnetic field together with the thermal FEM-equations. Each of them can be extended with an algebraic set of circuit equations. These include coupling terms as well, e.g. in resistances.

**The discretisations,** on which the computations for the single field problem are performed, do not have to be identical. Sometimes, only a sub-mesh has a physical meaning: e.g. air carrying a magnetic leakage flux is replaced by a convection constraint in the thermal model; the solid parts can be identical. Even the mesh in areas with more than one continuous degree of freedom can be discretised with different overlapping geometrical meshes and/or element types. Therefore, mesh transition operations have to be defined. The groups of algebraic equations can be solved with a strong coupled or with a cascade coupled strategy.

**The first approach consists of the generation** of a large system of non-linear equations with both types of FEM-equations, associated with the coupling terms. The mesh transitions and heat source terms have to be written as linearised algebraic functions. This large linear system may have unfavourable numerical properties due to the different nature of the underlying physical equations, resulting in a difficult to solve problem.

**The second method (Fig. 10.33)** defines an iterative process in which both sets of equations are solved sequentially. The mesh transitions and heat source calculations are necessary intermediate steps and do not necessarily have to be linearised. This approach can be interpreted as a "decomposition" of the different fields.

**Fig. 10.33. General flow chart of the cascade algorithm.**

**The following groups distinguish the operations necessary for the iterations:**

**• Updating of material properties:** The various temperature dependent material parameters are updated whenever the algorithm requires it. Therefore, various characteristics must be implemented. These involve electrical conductivities, thermal conductivities, permanent magnet material properties, thermal conductivities, loss coefficients and characteristics.

**• External process handling:** Calls to execute external FEM-solvers and mesh generators.

**• Iteration control:** Process commands to control the flow of the iterative loop and evaluation of stopping criteria. These criteria can be based on absolute or relative residuals or solution differences between two consecutive weighted solutions of the total problem or a sub-problem respectively.

**Fig. 10.34. Projection of element related values.**

**• Data transition commands:** When different meshes are used, mechanisms are necessary to project the field variables onto another mesh. Basically, this is a per-node interpolation of the solution.

**The projection of the values associated** with the element’s surface or volume, e.g. a calculated loss density, is not straightforward. If the meshes do not differ very much, the position of the centre of gravity of the element to be filled in, can be located in the other mesh (black dot Fig. 10.34). The corresponding element-related value can then be copied. If the meshes differ much or if a higher order of accuracy is desired, an average can be generated by means of a numerical integration using Gauss points (additional white dots Fig. 10.34).

**It can be stated that a large difference** between the meshes is not advantageous. It would mean that the related physical domains would not be calculated with a corresponding accuracy. However, mesh differences can arise due to local mesh quality reasons.

**• Adaptive relaxation of the convergence process:** In order to prevent the non-linear iteration process from divergence and to accelerate the convergence, an appropriate relaxation method must be applied. The damping factor can be predefined according to a certain function of the iteration number or adaptively, based on a minimisation of the total or partial (weighted) residual vector.

**• Heat sources calculation:** For the area covered by every meaningful element, a heat source density can be calculated based on the electromagnetic solution (Table 10.2).

**Table 10.2. Overview of main heat sources in electromagnetic problems.**

## Three-phase high voltage power cable

**In this section**, the results of a three-phase power cable simulation with respect to the coupled magnetic/electrostatic/thermal field problem are shown.

**Three-phase power cables exist in many variations and types, differing in conductor shape, material choice, conductor arrangement etc. They consist mainly of the following parts, in which several of the previously mentioned loss mechanisms can be found:**

**• conductor,** usually made of copper, suffering from joule losses caused by the high current

**• insulation layers and filling materials,** loaded with electric fields and therefore subject to dielectric losses

**Fig. 10.35. Geometry of the three-phase high voltage power cable.**

**• grounded lead sheath around** the primary isolation, shielding the electric field; due to its relative low conductance, internal eddy currents can develop

**• mechanical protections (armour),** sometimes made of magnetic steel and therefore subject to hysteresis and eddy current losses

**The presence of both electrically related** and magnetically related heat sources leads to a combined model consisting of three field types. Electric, magnetic and thermal fields that have to be calculated over a complete or a partial cross-section of the cable and its surrounding.

**The electrical field,** described by the scalar potential V, is only of interest in the isolation part loaded with an electrical field. Therefore, only a mesh covering this region is required to solve the electrostatic field equations.

The time-harmonic magnetic field is calculated on a larger mesh, since it is only partly shielded by the mechanical protection and thus a leakage field can exist outside the model. This leakage flux is considered by the region surrounding the cable geometry; the far field is modelled by a Kelvin transformed mesh. The losses consist of joule losses in the conducting regions, such as the lead, steel and copper and possible iron losses inside the steel.

**The thermal field is represented** by the temperature potential distribution T. The static thermal field region consists of the cable with the surrounding soil in which it is buried. From a certain distance, the ground is modelled by a Kelvin transformation and therefore assumed to be infinitely deep. It is assumed that the ground surface is cooled by convection.

**Fig. 10.36. Problem meshes used to compute the different physical fields.**

**The element-wise over the previous meshes** (electro-, magneto-static) calculated losses are projected onto the thermal mesh. The extracted temperatures are used to update the material properties in the other fields. Basically the dielectric loss factor and the thermal conductivity depend on temperature, but this is an effect of minor importance in this example problem. The largest parameter changes are encountered in the conductivity of the copper.

**To solve the entire problem,** a three-domain mesh must be constructed (Fig. 10.36). The mesh for the electric field contains 9943 first order triangular elements, the magnetic field mesh 15378 and the thermal field mesh 15560 elements.

The results of this threefold-coupled problem are collected in the following figures.

**Fig. 10.37. Entire solution of the threefold coupled problem.**

**The temperature in the centre** of a conductor amounts 84°C, which corresponds to reported measurements (Van Dommelen & Germay "6). Table 10.3 shows the computed values of the heat sources and their location within the models.

**Table 10.3. Overview of main heat sources in electromagnetic problems.**

## Coupled simulation for electrical machines

**An efficient simulation algorithm** for the coupled magnetic-thermal field of electric machines is proposed. This algorithm uses a combined FEM-circuit approach in both magnetic and thermal field regions. FEM calculations are performed to compute magnetic and thermal phenomena in a xy 2D Cartesian cross-section, whereas a circuit approach is applied to consider thermal phenomena in the z direction (Fig. 10.38).

**Fig. 10.38. Principle of the combined thermal FEM-circuit-model; the resistances in the model represent a part of the symmetry of the axial thermal phenomena.(The cross-links to the next slot are not shown.)**

The computations in the magnetic and the thermal domain result in four sets of equations: the magnetic FEM-equations, electrical circuit equations, the thermal FEM-equations and thermal circuit equations. The entire approach to solve the entire problem is a computationally efficient method. The simulation of a 15 kW TEFC (totally enclosed and fan-cooled) four-pole induction machine demonstrates the coupled approach.

**In general,** the heat transport perpendicular to the axis of cylindrical electrical machines can be modelled with a high accuracy. Less is known about the effects influencing the heat flow in the axial direction. For example, the flow of the cooling fluids in the end-regions is very complex. Therefore, 3D-calculations should be coupled to a ‘computational fluid dynamics’ simulation, asking for huge computational efforts.

The circuit approach offers the advantages that the heat transfer in the axial direction is modelled in a straightforward way, using thermal resistances. However, the determination of their values is troublesome and asks for experience and measurements.

**Different methods to solve the overall problem are possible:**

**• A cascade iteration algorithm in which** the systems are solved in successive steps. This results in a poor convergence. On the other hand loss calculation algorithms can be introduced in a simple way.

**• The solution of a numerically** strong coupled large system of equations in which the four sub-sets together with the coupling terms are assembled. This may result in a very ill-conditioned system of equations. Therefore, particular equation solvers are recommended to obtain a fast convergence.

**• An intermediate approach is possible** by placing the FEM-equations together with their corresponding circuit equations in one system of equations describing a single physical domain. These ‘physical’ sub-systems are then coupled in a cascade-like iteration. This method has a moderate rate of convergence, but requires less memory.

**The time-harmonic magnetic model is coupled with** circuit equations describing the effect of the end-windings, the bars outside the rotor and the end-rings as already introduced. The thermal model consists of a FEM component, including equivalent convection coefficients to account for ribs and equivalent thermal conductivities for areas such as the air gap and the slots. This thermal FEM model is extended by a circuit approach. The circuit equations represent the thermal paths connecting shaft, yokes, slots and the frame through thermal resistances and represent the internal end-region in air, the end-windings, end-rings, bearings and end-caps of the machine. The resulting temperature distribution of this coupled approach is plotted in Fig. 10.39. The isothermal lines in the shaft are caused by heat flowing through the shaft, a path that is described by a thermal network circuit equation.

**Fig. 10.39. Isothermal lines of the temperature distribution.**

### Modelling of thermal contact resistances

In the electromagnetic-thermal coupled modelling of electromagnetic devices, the thermal simulation poses some extra difficulties, not present in the electromagnetic field.

**Examples are the mixed** boundary conditions as found in convection (linear coefficient) or radiation (non-linear coefficient). These require extra adjustments of the finite-element matrix and right-hand-side.

**A second particularity is the presence** of thermal contact resistances, although it can be noted that electrical contact resistances exist as well when electrical current in the plane of the simulation is modelled. They are a particular problem encountered in thermal FEM models of the electromagnetic devices. Examples for this thermal contact interface are:

**• contact frame-stator yoke**

**• contact stator winding-slot**

**• contact rotor winding/bar-slot**

**• contact rotor yoke-shaft**

**• glue layer between a permanent magnet and the yoke.**

**Several approaches are known** to model such contact resistances. Corrections to the thermal conductivities of the conductor materials can be applied to consider this effect. This approach causes an error in the internal temperature distribution of the conducting region, but the average temperature is calculated accurately. The temperature distribution is used to update material data of the related sub-problems. Despite the extra computational costs of averaging, this approach has the advantage that there is an obvious geometrical relation between the elements of regions in the different sub-problems. This is advantageous for projection methods.

**An extra equivalent contact layer of elements can be inserted,** filled with an equivalent contact material. In order to obtain an acceptable aspect ratio of the finite elements, a sufficient number of elements have to be generated inside the contact layer and the adjacent regions. Due to the high number of slots in an electrical machine and the small size of the contact layer, this yields a significant growth of the numerical model and hence computation time. The number of elements can be reduced if the layer is enlarged, but this reduces the size of the conductor and the tooth, and so no clear geometrical relationship exists between the thermal and the magnetic sub-problem. If the same model, with the reduced conductor, were to be used for the magnetic sub-problem, an error would be made in the leakage flux and the joule heat calculation.

**Another possibility is to duplicate the** nodes lying on the edges marking the border between the conducting region and the iron. Extra terms in the equations modelling the contact are inserted to define the thermal relationship of the nodes. The standard meshing algorithms have to be adapted to generate the extra nodes. Moreover, these extra equations change the numerical properties of the matrix system to be solved.

**The third proposed approach will be outlined.** Here, the mesh generator must be able to double the nodes lying on the contact resistance (Fig. 10.40). This must be performed in a consistent way to prevent any ‘crossing’ of the edges of the contact resistance. Standard mesh generators need to be extended with appropriate search and data management structures.

**Fig. 10.40. Two elements at the thermal contact resistance interface.**

**Within the FEM code**, the standard element matrices are generated (the right-hand-side remains unaltered). The following extra matrix, suited for first order elements, is added to the system. The exact values of the coefficients are obtained using the Galerkin approach.

withthe contact resistance coefficient andthe length of the adjacent element edges.

A computed example using this approach is shown in Fig. 10.41. A partial thermal model of an induction machine is shown. The example consists of one stator winding slot and a rotor bar. The above relation models the contacts between the slots and the iron core. This leads to a set of isothermal lines appearing to be discontinuous. This is not true since many isothermal lines lie in the temperature jump inside the thermal contact resistance.

**Fig. 10.41. Computed temperature distribution using the thermal contact resistances in a slot arrangement of an electrical machine example.**