New Thermal-Conductivity Constitutive Matrix in Fourier’s Law for Heat Transfer Using the Cell Method

Featured Application

The amount of heat transferred by conduction is given by Fourier’s law. For the study of these phenomena, the application of computational techniques that allow the design of machines and devices used in engineering becomes crucial.

Abstract

In this paper, a new constitutive matrix $\left[M_{\tau }\right]$ for thermal conduction, for tetrahedral meshes, in a steady state thermal regime is developed through a new algebraic methodology, using the Cell Method as a computational method, which is included in the finite formulation. The constitutive matrix defines the behavior of solids when they are under a thermal potential. The results are compared with those obtained for the same problem by means of the constitutive matrix $\left[M_{\lambda }\right]$ developed previously, taking in both cases with a 2D axisymmetric model as reference, calculated with the finite element method. The errors obtained with the new matrix $\left[M_{\tau }\right]$ are of the order of 0.0025%, much lower than those obtained with the matrix $\left[M_{\lambda }\right]$.

1. Introduction

Solid metals have a high thermal conductivity. The transmission of heat by conduction is attributed to an exchange of energy between adjacent molecules and electrons in the conductive medium, without the macroscopic transfer of matter and without a visible displacement of particles.

The amount of heat transferred by conduction is given by Fourier’s law. This law states that the rate of heat conduction through a body per unit cross section is proportional to the temperature gradient that exists in the body.

For the study of these phenomena, the application of analytical methods and computational techniques that allow for the design of machines and devices used in engineering becomes very important.

The analytical methods proposed in this article are of an algebraic type. They have been applied to the study of heat transmission by conduction in a bimetallic tube. This case can be used for the study of pipes, or, in our particular case, for the future studies of the stator or rotor of an electric machine [1,2,3,4], as it can be seen in Figure 1. The most common heat sources in an electrical machine are the electromagnetic phenomena in the cores (approximated here as metal tubes), heated by the Joule effect in the conductors [5,6,7] and the mechanical friction in the moving parts [8,9].

In the present work, we propose the finite formulation (FF) [10,11,12,13,14,15], as well as the Cell Method (CM), [16,17] as an associated numerical method to analyze the numerical models proposed. In this methodology, we work with the global magnitudes associated with space-oriented elements such as volumes, surfaces, lines, and points of the discretized space, as well as to temporal elements, instead of the field magnitudes associated with independent variables with spatial and temporal coordinates [16,17,18,19,20].

In addition, the equations of constitutive type (equations of the medium) are clearly differentiated from the topological type (equations of balance) [10]. In FF, the physical laws that govern the thermal laws of heat transfer associated with electrical devices are expressed in their integral form. In this way, the final system of equations is posed directly, without the need to discretize the equivalent differential equations [20].

The thermal analysis of electrical devices using this methodology greatly facilitates the implementation of the boundary conditions and continuity when working with global magnitudes, and, furthermore, directly raises the system of equations without the need to discretize the differential equations.

Three methodologies have been formulated previously in order to obtain the constitutive thermal matrix. These methods are those by Tonti [16], Bullo [21,22], and the method proposed by Specogna [14], for problems of electrical conduction, which we have adapted to thermal conduction [10,11,12,13].

These three methods use the projections of edges and surfaces from dual to primal space. They also use local coordinate systems with subsequent transformations to global coordinates. We calculate the barycenters of the dual surfaces, and then obtain a weighted dual barycenter from those that were previously found.

Tonti [16] discusses a new 2D numerical method for the solution of temperature field equations using CM. The essence of the method is to directly provide a discrete formulation of field laws. It is proven that, for linear interpolation, the stiffness matrix obtained coincides with the one of the finite element method (FEM). For quadratic interpolation, however, the present stiffness matrix differs from that of FEM; moreover, it is asymmetric. It is shown that by using a parabolic interpolation, a convergence of the fourth order is obtained. This is greater than the one obtained with FEM, using the same interpolation.

Bullo [21] calculates, through CM, the 2D fields applied to a coupled computation of electric and thermal conduction using a linear interpolation of both the electric and temperature fields, and uses a quadratic interpolation for the thermal analysis approach. Bullo [22] uses CM for the solution of coupled problems of steady-state electric and transient thermal conductions in 3D regions. Dual barycentric cell complexes are used for both space and time domains, the latter inducing a Crank–Nicolson time integration scheme.

Specogna [14], by using a CM formulation for eddy currents, presents a geometric approach to construct approximations of the discrete magnetic and Ohm’s constitutive matrices. In the case of the Ohm’s matrix, he also shows how to make it symmetric. He compares the impact on the solution of the proposed Ohm’s matrices, and an iterative technique to obtain a consistent right-hand-side term in the final system is described.

González [12,13] developed a new constitutive matrix ($\left[M_{\lambda }\right]$) for thermal conduction in a transient thermal regime using CM. He demonstrates that this matrix is equivalent to the electrical conduction constitutive matrix in a steady state, and applies this constitutive matrix to the thermal analysis of asynchronous electric machines in a transient regime.

Monzón-Verona [12,13] analyzed the temperature distribution in a conductor disk in a transitory regime. The disk is in motion in a stationary magnetic field generated by a permanent magnet, and so the electric currents induced inside it generate heat. The system acts as a magnetic brake, and is analyzed using infrared sensor techniques. In addition, for the simulation and analysis of the magnetic brake, a new thermal convective matrix for the 3D Cell Method (CM) is proposed.

In the present article, a new constitutive thermal conduction matrix $\left[M_{\tau }\right]$ is obtained with better results than with $\left[M_{\lambda }\right]$.

The new thermal conduction constitutive matrix $\left[M_{\tau }\right]$ formulated with CM in 3D has been verified by contrasting the numerical results with those obtained by FEM. The difference between the CM results obtained and those obtained in FEM is less than 0.0025%.

The main advantage of the method proposed in this article is its simplicity. The constitutive matrices developed by previous methods presented complex calculations, while the new constitutive matrix depends exclusively on the coordinates of the vertices of the tetrahedra, which constitute the mesh.

This work has been divided into the following sections: Section 2 explains, in detail, the analytical methodology for obtaining the new thermal conduction constitutive matrix $\left[M_{\tau }\right]$ in CM, formulating the corresponding conductive term. Section 3 shows the results obtained, validating the previous formulation through the computational simulation with $\left[M_{\lambda }\right]$ and $\left[M_{\tau }\right]$. Finally, Section 4 presents the conclusions.

2. Calculation of the New Thermal Constitutive Matrix $\left[{\mathit{M}}_{\mathit{\tau }}\right]$

In the present section, we obtain the analytical formulation of the matrix $\left[M_{\tau }\right]$ that has been developed to improve the results obtained with $\left[M_{\lambda }\right]$. In a steady state, and without internal heat sources, the equation of energy balance without mass transfer [10,11,12,13] is in the CM,

$$\tilde{D}\left(-M_{\tau }GT\right) = G^t\left(M_{\tau }GT\right) = 0.$$

(1)

The domain is meshed by tetrahedral elements. In CM, the tetrahedron of Figure 2 will be taken as the reference cell, with its nodes, and edges and dual surfaces oriented internally and externally, respectively.

In CM, the discrete gradient operator for the reference tetrahedron is defined as follows:

$$\left[G\right] = \left[\begin{array}{cc}\begin{array}{cc}-1& 1\\ -1& 0\\ -1& 0\end{array}& \begin{array}{cc}0& 0\\ 1& 0\\ 0& 1\end{array}\\ \begin{array}{cc}0& -1\\ 0& 0\\ 0& -1\end{array}& \begin{array}{cc}1& 0\\ -1& 1\\ 0& 1\end{array}\end{array}\right].$$

(2)

It is assumed that the temperature is distributed within the tetrahedron, following the affine function of Cartesian spatial coordinates

$${\tau }_i\left(x, y, z\right) = g_x x+g_y y+g_z z+a,$$

(3)

where $a$ is an auxiliary constant introduced in order to later develop a square matrix. Then, the temperatures in the primary nodes of the reference tetrahedron are as follows:

$$\begin{array}{c}\begin{array}{c}{\tau }_0 = g_x x_0+g_y y_0+g_z z_0+a\\ {\tau }_1 = g_x x_1+g_y y_1+g_z z_1+a\end{array}\\ \begin{array}{c}{\tau }_2 = g_x x_2+g_y y_2+g_z z_2+a\\ {\tau }_3 = g_x x_3+g_y y_3+g_z z_3+a\end{array}\end{array},$$

(4)

where the Cartesian coordinates of the nodes $N = \left\{n_0,n_1,n_2,n_3\right\}$ are the following

$$\left[\begin{array}{c}\begin{array}{c}{n}_0\\ n_1\end{array}\\ \begin{array}{c}{n}_2\\ n_3\end{array}\end{array}\right] = \left[\begin{array}{c}\begin{array}{c}\begin{array}{ccc}{x}_0& y_0& z_0\end{array}\\ \begin{array}{ccc}{x}_1& y_1& z_1\end{array}\end{array}\\ \begin{array}{c}\begin{array}{ccc}{x}_2& y_2& z_2\end{array}\\ \begin{array}{ccc}{x}_3& y_3& z_3\end{array}\end{array}\end{array}\right].$$

(5)

The equation system shown in Equation (4), in matrix form, is as follows:

$$\left[\begin{array}{cc}\begin{array}{cc}\begin{array}{c}\begin{array}{c}{x}_0\\ x_1\end{array}\\ \begin{array}{c}{x}_2\\ x_3\end{array}\end{array}& \begin{array}{c}\begin{array}{c}{y}_0\\ y_1\end{array}\\ \begin{array}{c}{y}_2\\ y_3\end{array}\end{array}\end{array}& \begin{array}{cc}\begin{array}{c}\begin{array}{c}{z}_0\\ z_1\end{array}\\ \begin{array}{c}{z}_2\\ z_3\end{array}\end{array}& \begin{array}{c}\begin{array}{c}1\\ 1\end{array}\\ \begin{array}{c}1\\ 1\end{array}\end{array}\end{array}\end{array}\right]\left[\begin{array}{c}\begin{array}{c}{g}_x\\ g_y\end{array}\\ \begin{array}{c}{g}_z\\ a\end{array}\end{array}\right] = \left[\begin{array}{c}\begin{array}{c}{\tau }_0\\ {\tau }_1\end{array}\\ \begin{array}{c}{\tau }_2\\ {\tau }_3\end{array}\end{array}\right],$$

(6)

which in short form can be written as

$${\left[B\right]}_{4\times 4}{\left[G_a\right]}_{4\times 1} = {\left[\tau \right]}_{4\times 1},$$

(7)

and then,

$$\left[G_a\right] = {\left[B\right]}^{-1}\left[\tau \right].$$

(8)

The Cramer rule can be used to solve Equation (6). The determinant of the system is

$$\Delta = \left|\begin{array}{cc}\begin{array}{cc}\begin{array}{c}\begin{array}{c}{x}_0\\ x_1\end{array}\\ \begin{array}{c}{x}_2\\ x_3\end{array}\end{array}& \begin{array}{c}\begin{array}{c}{y}_0\\ y_1\end{array}\\ \begin{array}{c}{y}_2\\ y_3\end{array}\end{array}\end{array}& \begin{array}{cc}\begin{array}{c}\begin{array}{c}{z}_0\\ z_1\end{array}\\ \begin{array}{c}{z}_2\\ z_3\end{array}\end{array}& \begin{array}{c}\begin{array}{c}1\\ 1\end{array}\\ \begin{array}{c}1\\ 1\end{array}\end{array}\end{array}\end{array}\right|.$$

(9)

Then,

$$g_x = \frac{\left|\begin{array}{cc}\begin{array}{cc}\begin{array}{c}\begin{array}{c}{\tau }_0\\ {\tau }_1\end{array}\\ \begin{array}{c}{\tau }_2\\ {\tau }_3\end{array}\end{array}& \begin{array}{c}\begin{array}{c}{y}_0\\ y_1\end{array}\\ \begin{array}{c}{y}_2\\ y_3\end{array}\end{array}\end{array}& \begin{array}{cc}\begin{array}{c}\begin{array}{c}{z}_0\\ z_1\end{array}\\ \begin{array}{c}{z}_2\\ z_3\end{array}\end{array}& \begin{array}{c}\begin{array}{c}1\\ 1\end{array}\\ \begin{array}{c}1\\ 1\end{array}\end{array}\end{array}\end{array}\right|}{\Delta }, g_y = \frac{\left|\begin{array}{cc}\begin{array}{cc}\begin{array}{c}\begin{array}{c}{x}_0\\ x_1\end{array}\\ \begin{array}{c}{x}_2\\ x_3\end{array}\end{array}& \begin{array}{c}\begin{array}{c}{\tau }_0\\ {\tau }_1\end{array}\\ \begin{array}{c}{\tau }_2\\ {\tau }_3\end{array}\end{array}\end{array}& \begin{array}{cc}\begin{array}{c}\begin{array}{c}{z}_0\\ z_1\end{array}\\ \begin{array}{c}{z}_2\\ z_3\end{array}\end{array}& \begin{array}{c}\begin{array}{c}1\\ 1\end{array}\\ \begin{array}{c}1\\ 1\end{array}\end{array}\end{array}\end{array}\right|}{\Delta }, g_z = \frac{\left|\begin{array}{cc}\begin{array}{cc}\begin{array}{c}\begin{array}{c}{x}_0\\ x_1\end{array}\\ \begin{array}{c}{x}_2\\ x_3\end{array}\end{array}& \begin{array}{c}\begin{array}{c}{y}_0\\ y_1\end{array}\\ \begin{array}{c}{y}_2\\ y_3\end{array}\end{array}\end{array}& \begin{array}{cc}\begin{array}{c}\begin{array}{c}{\tau }_0\\ {\tau }_1\end{array}\\ \begin{array}{c}{\tau }_2\\ {\tau }_3\end{array}\end{array}& \begin{array}{c}\begin{array}{c}1\\ 1\end{array}\\ \begin{array}{c}1\\ 1\end{array}\end{array}\end{array}\end{array}\right|}{\Delta }.$$

(10)

We do not calculate $a$ because $\raisebox{1ex}{$\partial a$}\!\left/ \!\raisebox{-1ex}{$\partial x$}\right. = \raisebox{1ex}{$\partial a$}\!\left/ \!\raisebox{-1ex}{$\partial y $}\right.=\raisebox{1ex}{$\partial a$}\!\left/ \!\raisebox{-1ex}{$\partial z$}\right. = 0$.

2.1. Analytical Development of g_x, g_y, and g_z

To calculate g_x, we developed Equation (10) using adjuncts of the determinant

$$g_x = {\tau }_0\frac{\left[\begin{array}{ccc}{y}_1& z_1& 1\\ y_2& z_2& 1\\ y_3& z_3& 1\end{array}\right]}{\Delta }-{\tau }_1\frac{\left[\begin{array}{ccc}{y}_0& z_0& 1\\ y_2& z_2& 1\\ y_3& z_3& 1\end{array}\right]}{\Delta }+{\tau }_2\frac{\left[\begin{array}{ccc}{y}_0& z_0& 1\\ y_1& z_1& 1\\ y_3& z_3& 1\end{array}\right]}{\Delta }-{\tau }_3\frac{\left[\begin{array}{ccc}{y}_0& z_0& 1\\ y_1& z_1& 1\\ y_2& z_2& 1\end{array}\right]}{\Delta },$$

(11)

and we obtain the following expression

$$\begin{array}{cc}{g}_x = & \frac{{\tau }_0}{\Delta }\left[\left(y_2{z}_3-y_3{z}_2\right)-\left(y_1{z}_3-y_3{z}_1\right)+\left(y_1{z}_2-y_2{z}_1\right)\right]\\ & -\frac{{\tau }_1}{\Delta }\left[\left(y_2{z}_3-y_3{z}_2\right)-\left(y_0{z}_3-y_3{z}_0\right)+\left(y_0{z}_2-y_2{z}_0\right)\right]\\ & +\frac{{\tau }_2}{\Delta }\left[\left(y_1{z}_3-y_3{z}_1\right)-\left(y_0{z}_3-y_3{z}_0\right)+\left(y_0{z}_1-y_1{z}_0\right)\right]\\ & -\frac{{\tau }_3}{\Delta }\left[\left(y_1{z}_2-y_2{z}_1\right)-\left(y_0{z}_2-y_2{z}_0\right)+\left(y_0{z}_1-y_1{z}_0\right)\right]\end{array}.$$

(12)

In the same way, we obtain g_y and y g_z

$$\begin{array}{cc}{g}_y = & -\frac{{\tau }_0}{\Delta }\left[\left(x_2{z}_3-x_3{z}_2\right)-\left(x_1{z}_3-x_3{z}_1\right)+\left(x_1{z}_2-x_2{z}_1\right)\right]\\ & +\frac{{\tau }_1}{\Delta }\left[\left(x_2{z}_3-x_3{z}_2\right)-\left(x_0{z}_3-x_3{z}_0\right)+\left(x_0{z}_2-x_2{z}_0\right)\right]\\ & -\frac{{\tau }_2}{\Delta }\left[\left(x_1{z}_3-x_3{z}_1\right)-\left(x_0{z}_3-x_3{z}_0\right)+\left(x_0{z}_1-x_1{z}_0\right)\right]\\ & +\frac{{\tau }_3}{\Delta }\left[\left(x_1{z}_2-x_2{z}_1\right)-\left(x_0{z}_2-x_2{z}_0\right)+\left(x_0{z}_1-x_1{z}_0\right)\right]\end{array},$$

(13)

$$\begin{array}{cc}{g}_z = & \frac{{\tau }_0}{\Delta }\left[\left(x_2{z}_3-x_3{y}_2\right)-\left(x_1{y}_3-x_3{y}_1\right)+\left(x_1{y}_2-x_2{y}_1\right)\right]\\ & -\frac{{\tau }_1}{\Delta }\left[\left(x_2{y}_3-x_3{y}_2\right)-\left(x_0{y}_3-x_3{y}_0\right)+\left(x_0{y}_2-x_2{y}_0\right)\right]\\ & +\frac{{\tau }_2}{\Delta }\left[\left(x_1{y}_3-x_3{y}_1\right)-\left(x_0{y}_3-x_3{y}_0\right)+\left(x_0{y}_1-x_1{y}_0\right)\right]\\ & -\frac{{\tau }_3}{\Delta }\left[\left(x_1{y}_2-x_2{y}_1\right)-\left(x_0{y}_2-x_2{y}_0\right)+\left(x_0{y}_1-x_1{y}_0\right)\right]\end{array}.$$

(14)

We know that

$$grad \tau \left(x, y, z\right) = \frac{\partial \tau }{\partial x}\overrightarrow{i}+\frac{\partial \tau }{\partial y}\overrightarrow{j}+\frac{\partial \tau }{\partial z}\overrightarrow{k}.$$

(15)

Taking into account Equation (2), then,

$$\frac{\partial \tau }{\partial x} = g_x, \frac{\partial \tau }{\partial y} = g_y, \frac{\partial \tau }{\partial z} = g_z,$$

(16)

and, hence,

$$grad \tau \left(x, y, z\right) = g_x \overrightarrow{i}+g_{y }\overrightarrow{j}+g_z \overrightarrow{k}.$$

(17)

2.2. Building the Matrix $\left[A_{\tau }\right]$

The Fourier heat transmission equation has been established, using the CM, as follows

$${\left[Q^a\right]}_{6\times 1} = {\left[M_{\tau }\right]}_{6\times 6}{\left[G\right]}_{6\times 4}{\left[\tau \right]}_{4\times 1},$$

(18)

where $\left[Q^a\right]$ is the heat flow transmitted and $\left[M_{\tau }\right]$ is the constitutive matrix of thermal transmission we propose in this article. Suppose there is a matrix $\left[A_{\tau }\right]$, such that

$${\left[M_{\tau }\right]}_{6\times 6} = {\left[\tilde{S}\right]}_{6\times 3}{\left[A_{\tau }\right]}_{3\times 6},$$

(19)

where $\left[\tilde{S}\right]$ is the dual faces matrix of the cell

$$\left[\tilde{S}\right] = \left[\begin{array}{ccc}\begin{array}{c}\begin{array}{c}{S}_{0x}\\ S_{1x}\\ S_{2x}\end{array}\\ \begin{array}{c}{S}_{3x}\\ S_{4x}\\ S_{5x}\end{array}\end{array}& \begin{array}{c}\begin{array}{c}{S}_{0y}\\ S_{1y}\\ S_{2y}\end{array}\\ \begin{array}{c}{S}_{3y}\\ S_{4y}\\ S_{5y}\end{array}\end{array}& \begin{array}{c}\begin{array}{c}{S}_{0z}\\ S_{1z}\\ S_{2z}\end{array}\\ \begin{array}{c}{S}_{3z}\\ S_{4z}\\ S_{5z}\end{array}\end{array}\end{array}\right] .$$

(20)

The heat flow transmitted $\left[q\right]$ is

$${\left[q\right]}_{3\times 1} = {\left[A_{\tau }\right]}_{3\times 6}{\left[G\right]}_{6\times 4}{\left[\tau \right]}_{4\times 1}.$$

(21)

We define an unknown vector $\left[X\right]$, such as

$${\left[X\right]}_{6\times 1} = {\left[G\right]}_{6\times 4}{\left[\tau \right]}_{4\times 1}.$$

(22)

The density heat vector $\left[\overrightarrow{q}\right]$ is

$${\left[\overrightarrow{q}\right]}_{3\times 1} = -\lambda \left[\begin{array}{c}{g}_x\\ g_y\\ g_z\end{array}\right] = -\lambda {\left[grad\left(\tau \right)\right]}_{3\times 1},$$

(23)

where $\lambda$ is the thermal conductivity coefficient, and, hence, the heat flow $\left[Q^b\right]$ is

$${\left[Q^b\right]}_{6\times 1} = {\left[\tilde{S}\right]}_{6\times 3}{\left[\overrightarrow{q}\right]}_{3\times 1}.$$

(24)

Then, what is written in Equation (23) is equal to Equation (18), and therefore $\left[Q^b\right] = \left[Q^a\right]$

$$-\lambda {\left[\tilde{S}\right]}_{6\times 3}{\left[grad\left(\tau \right)\right]}_{3\times 1} = {\left[M_{\tau }\right]}_{6\times 6}{\left[G\right]}_{6\times 4}{\left[\tau \right]}_{4\times 1}.$$

(25)

Replacing Equation (18) and Equation (21) in the second member of Equation (25), we obtain

$$-\lambda {\left[\tilde{S}\right]}_{6\times 3}{\left[grad\left(\tau \right)\right]}_{3\times 1} = {\left[\tilde{S}\right]}_{6\times 3}{\left[A_{\tau }\right]}_{3\times 6}{\left[X\right]}_{6\times 1}.$$

(26)

Simplifying the dual faces matrix $\left[\tilde{S}\right]$

$$-\lambda {\left[grad\left(\tau \right)\right]}_{3\times 1} = {\left[A_{\tau }\right]}_{3\times 6}{\left[X\right]}_{6\times 1}.$$

(27)

Replacing $\left[X\right]$ with the values of $\left[G\right]$ and $\left[\tau \right]$, we obtain

$${\left[X\right]}_{6\times 1} = {\left[G\right]}_{6\times 4}{\left[\tau \right]}_{4\times 1} = {\left[\begin{array}{cc}\begin{array}{cc}-1& 1\\ -1& 0\\ -1& 0\end{array}& \begin{array}{cc}0& 0\\ 1& 0\\ 0& 1\end{array}\\ \begin{array}{cc}0& -1\\ 0& 0\\ 0& -1\end{array}& \begin{array}{cc} 1& 0\\ -1& 1\\ 0& 1\end{array}\end{array}\right]}_{6\times 4} {\left[\begin{array}{c}\begin{array}{c}{\tau }_0\\ {\tau }_1\end{array}\\ \begin{array}{c}{\tau }_2\\ {\tau }_3\end{array}\end{array}\right]}_{4\times 1} = {\left[\begin{array}{c}\begin{array}{c}-{\tau }_0+{\tau }_1\\ -{\tau }_0+{\tau }_2\\ -{\tau }_0+{\tau }_3\end{array}\\ \begin{array}{c}-{\tau }_1+{\tau }_2\\ -{\tau }_2+{\tau }_3\\ -{\tau }_1+{\tau }_3\end{array}\end{array}\right]}_{6\times 1}$$

(28)

and, therefore,

$${\left[\begin{array}{c}\begin{array}{c}{x}_0\\ x_1\\ x_2\end{array}\\ \begin{array}{c}{x}_3\\ x_4\\ x_5\end{array}\end{array}\right]}_{6\times 1} = {\left[\begin{array}{c}\begin{array}{c}-{\tau }_0+{\tau }_1\\ -{\tau }_0+{\tau }_2\\ -{\tau }_0+{\tau }_3\end{array}\\ \begin{array}{c}-{\tau }_1+{\tau }_2\\ -{\tau }_2+{\tau }_3\\ -{\tau }_1+{\tau }_3\end{array}\end{array}\right]}_{6\times 1}.$$

(29)

Then, observing Equation (27), replacing the value of ${\left[grad\left(\tau \right)\right]}_{3\times 1}$, then,

$$-\lambda {\left[\begin{array}{c}{g}_x\\ g_y\\ g_z\end{array}\right]}_{3\times 1} = {\left[A_{\tau }\right]}_{3\times 6}{\left[\begin{array}{c}\begin{array}{c}{x}_0\\ x_1\\ x_2\end{array}\\ \begin{array}{c}{x}_3\\ x_4\\ x_5\end{array}\end{array}\right]}_{6\times 1},$$

(30)

that is to say

$$-\lambda {\left[\begin{array}{c}{g}_x\\ g_y\\ g_z\end{array}\right]}_{3\times 1} = {\left[\begin{array}{cc}\begin{array}{ccc}\begin{array}{c}{A}_{00}\\ A_{10}\\ A_{20}\end{array}& \begin{array}{c}{A}_{01}\\ A_{11}\\ A_{21}\end{array}& \begin{array}{c}{A}_{02}\\ A_{12}\\ A_{22}\end{array}\end{array}& \begin{array}{ccc}\begin{array}{c}{A}_{03}\\ A_{13}\\ A_{23}\end{array}& \begin{array}{c}{A}_{04}\\ A_{14}\\ A_{24}\end{array}& \begin{array}{c}{A}_{05}\\ A_{15}\\ A_{25}\end{array}\end{array}\end{array}\right]}_{3\times 6}{\left[\begin{array}{c}\begin{array}{c}{x}_0\\ x_1\\ x_2\end{array}\\ \begin{array}{c}{x}_3\\ x_4\\ x_5\end{array}\end{array}\right]}_{6\times 1},$$

(31)

and, therefore, developing Equation (31), we obtain the following expression

$$-\lambda \left[\begin{array}{c}{g}_x\\ g_y\\ g_z\end{array}\right] = \left[\begin{array}{c}{A}_{00} x_0+A_{01} x_1+A_{02} x_2+A_{03} x_3+A_{04} x_4+A_{05} x_5\\ A_{10} x_0+A_{11} x_1+A_{12} x_2+A_{13} x_3+A_{14} x_4+A_{15} x_5\\ A_{20} x_0+A_{21} x_1+A_{22} x_2+A_{23} x_3+A_{24} x_4+A_{25} x_5\end{array}\right].$$

(32)

The value of $g_x$ is included in the following equation

$$-\lambda g_x = A_{00} x_0+A_{01} x_1+A_{02} x_2+A_{03} x_3+A_{04} x_4+A_{05} x_5.$$

(33)

Replacing the values $x_i$ developed in Equation (29), we obtain

$$\begin{gathered} -\lambda g_x = A_{00} \left(-{\tau }_0+{\tau }_1\right)+A_{01} \left(-{\tau }_0+{\tau }_2\right)+A_{02} \left(-{\tau }_0+{\tau }_3\right)+ \\ +A_{03} \left(-{\tau }_1+{\tau }_2\right)+A_{04} \left(-{\tau }_2+{\tau }_3\right)+A_{05} \left(-{\tau }_1+{\tau }_3\right). \end{gathered}$$

(34)

If we group the terms affected by the same temperature value, then

$$\begin{array}{cc}\begin{array}{c}-\lambda g_x = \\ \begin{array}{c}\\ \begin{array}{c}\\ \end{array}\end{array}\end{array}& \begin{array}{c}{\tau }_0\left(-A_{00}-A_{01}-A_{02}\right)+\\ \begin{array}{c}+{\tau }_1\left(A_{00}-A_{03}-A_{05}\right)+\\ \begin{array}{c}+{\tau }_2\left(A_{01}+A_{03}-A_{04}\right)+\\ +{\tau }_3\left(A_{02}+A_{04}+A_{05}\right) , \end{array}\end{array}\end{array}\end{array}$$

(35)

and replacing the value of $g_x$ calculated in Equation (12), then, by matching what was obtained in Equation (35), we get

$$\begin{array}{cc}-\frac{\lambda {\tau }_0}{\Delta }\left[\left(y_2{z}_3-y_3{z}_2\right)-\left(y_1{z}_3-y_3{z}_1\right)+\left(y_1{z}_2-y_2{z}_1\right)\right]& \\ +\frac{\lambda {\tau }_1}{\Delta }\left[\left(y_2{z}_3-y_3{z}_2\right)-\left(y_0{z}_3-y_3{z}_0\right)+\left(y_0{z}_2-y_2{z}_0\right)\right]& \\ -\frac{\lambda {\tau }_2}{\Delta }\left[\left(y_1{z}_3-y_3{z}_1\right)-\left(y_0{z}_3-y_3{z}_0\right)+\left(y_0{z}_1-y_1{z}_0\right)\right]& \\ +\frac{\lambda {\tau }_3}{\Delta }\left[\left(y_1{z}_2-y_2{z}_1\right)-\left(y_0{z}_2-y_2{z}_0\right)+\left(y_0{z}_1-y_1{z}_0\right)\right] = & \\ & {\tau }_0\left(-A_{00}-A_{01}-A_{02}\right)\\ & +{\tau }_1\left(A_{00}-A_{03}-A_{05}\right)\\ & +{\tau }_2\left(A_{01}+A_{03}-A_{04}\right)\\ & +{\tau }_3\left(A_{02}+A_{04}+A_{05}\right)\end{array}.$$

(36)

Then, matching the terms that affect to ${\tau }_0$ in (36) in both sides, we obtain the following expression

$$-\frac{\lambda {\tau }_0}{\Delta }\left[\left(y_2{z}_3-y_3{z}_2\right)-\left(y_1{z}_3-y_3{z}_1\right)+\left(y_1{z}_2-y_2{z}_1\right)\right] = {\tau }_0\left(-A_{00}-A_{01}-A_{02}\right).$$

(37)

Comparing the terms in Equation (37)

$$\begin{array}{ccc} -\frac{\lambda }{\Delta }\left(y_2{z}_3-y_3{z}_2\right)& =& -A_{00}\\ \frac{\lambda }{\Delta }\left(y_1{z}_3-y_3{z}_1\right)& =& -A_{01}\\ -\frac{\lambda }{\Delta }\left(y_1{z}_2-y_2{z}_1\right)& =& -A_{02}\end{array} \},$$

(38)

we get

$$\begin{array}{c}{A}_{00} = \frac{\lambda }{\Delta }\left(y_2{z}_3-y_3{z}_2\right)\\ A_{01} = \frac{\lambda }{\Delta }\left(y_3{z}_1-y_1{z}_3\right)\\ A_{02} = \frac{\lambda }{\Delta }\left(y_1{z}_2-y_2{z}_1\right)\end{array}.$$

(39)

Matching the terms that affect to ${\tau }_1$ in Equation (36)

$$\frac{\lambda {\tau }_1}{\Delta }\left[\left(y_2{z}_3-y_3{z}_2\right)-\left(y_0{z}_3-y_3{z}_0\right)+\left(y_0{z}_2-y_2{z}_0\right)\right] = {\tau }_1\left(A_{00}-A_{03}-A_{05}\right),$$

(40)

we obtain

$$\begin{array}{c}{A}_{00} = \frac{\lambda }{\Delta }\left(y_2{z}_3-y_3{z}_2\right)\\ A_{03} = \frac{\lambda }{\Delta }\left(y_0{z}_3-y_3{z}_0\right)\\ A_{05} = \frac{\lambda }{\Delta }\left(y_2{z}_0-y_0{z}_2\right)\end{array}.$$

(41)

Matching the terms that affect to ${\tau }_2$ in Equation (36), then,

$$-\frac{\lambda {\tau }_2}{\Delta }\left[\left(y_1{z}_3-y_3{z}_1\right)-\left(y_0{z}_3-y_3{z}_0\right)+\left(y_0{z}_1-y_1{z}_0\right)\right] = +{\tau }_2\left(A_{01}+A_{03}-A_{04}\right),$$

(42)

and, therefore, we get

$$\begin{array}{c}{A}_{01}=\frac{\lambda }{\Delta }\left(y_3{z}_1-y_1{z}_3\right)\\ A_{03}=\frac{\lambda }{\Delta }\left(y_0{z}_3-y_3{z}_0\right)\\ A_{04}=\frac{\lambda }{\Delta }\left(y_0{z}_1-y_1{z}_0\right)\end{array}.$$

(43)

In the same way, matching the terms that affect to ${\tau }_3$ in Equation (36), then,

$$\frac{\lambda {\tau }_3}{\Delta }\left[\left(y_1{z}_2-y_2{z}_1\right)-\left(y_0{z}_2-y_2{z}_0\right)+\left(y_0{z}_1-y_1{z}_0\right)\right]={\tau }_3\left(A_{02}+A_{04}+A_{05}\right),$$

(44)

and, hence, we get

$$\begin{array}{c}{A}_{02}=\frac{\lambda }{\Delta }\left(y_1{z}_2-y_2{z}_1\right)\\ A_{05}=\frac{\lambda }{\Delta }\left(y_2{z}_0-y_0{z}_2\right)\\ A_{04}=\frac{\lambda }{\Delta }\left(y_0{z}_1-y_1{z}_0\right)\end{array}.$$

(45)

Following the same procedure, we obtain the terms of $g_y$

$$\begin{array}{c}{A}_{10} = \frac{\lambda }{\Delta }\left(x_3{z}_2-x_2{z}_3\right)\\ A_{11} = \frac{\lambda }{\Delta }\left(x_1{z}_3-x_3{z}_1\right)\\ A_{12} = \frac{\lambda }{\Delta }\left(x_2{z}_1-x_1{z}_2\right)\end{array}\hspace{0.17em}\begin{array}{c}{A}_{13} = \frac{\lambda }{\Delta }\left(x_3{z}_0-x_0{z}_3\right)\\ A_{14} = \frac{\lambda }{\Delta }\left(x_1{z}_0-x_0{z}_1\right)\\ A_{15} = \frac{\lambda }{\Delta }\left(x_0{z}_2-x_2{z}_0\right)\end{array},$$

(46)

and of $g_z$

$$\begin{array}{c}{A}_{20}=\frac{\lambda }{\Delta }\left(x_2{z}_3-x_3{y}_2\right)\\ A_{21}=\frac{\lambda }{\Delta }\left(x_3{y}_1-x_1{y}_3\right)\\ A_{22}=\frac{\lambda }{\Delta }\left(x_1{y}_2-x_2{y}_1\right)\end{array}\hspace{0.17em}\begin{array}{c}{A}_{23}=\frac{\lambda }{\Delta }\left(x_0{y}_3-x_3{y}_0\right)\\ A_{24}=\frac{\lambda }{\Delta }\left(x_0{y}_1-x_1{y}_0\right)\\ A_{25}=\frac{\lambda }{\Delta }\left(x_2{y}_0-x_0{y}_2\right)\end{array}.$$

(47)

2.3. New $\left[M_{\tau }\right]$ Matrix

In the previous section, the terms for A_ij have been obtained. Now we can build the matrix $\left[A_{\tau }\right]$,

$$[A_{\tau }] = \frac{\lambda }{\Delta }\left[\begin{array}{cccccc}\left(y_3{z}_2-y_2{z}_3\right)& \left(y_1{z}_3-y_3{z}_1\right)& \left(y_2{z}_1-y_1{z}_2\right)& \left(y_3{z}_0-y_0{z}_3\right)& \left(y_1{z}_0-y_0{z}_1\right)& \left(y_0{z}_2-y_2{z}_0\right)\\ \left(x_2{z}_3-x_3{z}_2\right)& \left(x_3{z}_1-x_1{z}_3\right)& \left(x_1{z}_2-x_2{z}_1\right)& \left(x_0{z}_3-x_3{z}_0\right)& \left(x_0{z}_1-x_1{z}_0\right)& \left(x_2{z}_0-x_0{z}_2\right)\\ \left(x_3{y}_2-x_2{y}_3\right)& \left(x_1{y}_3-x_3{y}_1\right)& \left(x_2{y}_1-x_1{y}_2\right)& \left(x_3{y}_0-x_0{y}_3\right)& \left(x_1{y}_0-x_0{y}_1\right)& \left(x_0{y}_2-x_2{y}_0\right)\end{array}\right].$$

(48)

Therefore, the new thermal conductivity constitutive matrix is

$${\left[M_{\tau }\right]}_{6\times 6} = {\left[\tilde{S}\right]}_{6\times 3}{\left[A_{\tau }\right]}_{3\times 6}.$$

(49)

3. Results and Validation

Verification consists of checking that the procedure implemented is conceptually correct, and validation consists of checking that the data obtained from the numerical simulations coincide with an objective reality [23,24]. The verification was carried out in Section 2. To check the validation of the new matrix, three types of numerical simulations were designed and executed, which will be explained below. On the one hand, a highly accurate numerical simulation was carried out using FEM, and then two numerical simulations were carried out to evaluate the accuracy of $\left[M_{\lambda }\right]$ and $\left[M_{\tau }\right]$, using CM.

3.1. Validation of FEM

Our reference to validate the results obtained in Section 3.2 and Section 3.3 was FEM. Specifically, we used the FEMM program (Finite Element Method Magnetics program) [25], and then we applied the statistics to verify FEMM´s validity.

We added this section to verify FEMM through the analytical solution of a simple problem for a single thermal conductivity. Then, FEMM was applied as a reference tool to verify a more complex problem of two thermal conductivities, and we compared the FEM with the results obtained by means of the two thermal conductivity matrices, $\left[M_{\lambda }\right]$ and $\left[M_{\tau }\right]$, analyzed with the CM.

The problem consists of analyzing the distribution of temperatures in a tube with a single thermal conductivity and axial symmetry. In Figure 3, a section of the tube with dimensions 0.5 × 2.0 m can be observed, associated to a cylinder 2 m high and 2 m in diameter. The cylinder wall had a thermal conduction coefficient of 1 Wm⁻¹K⁻¹.

This problem has the following analytical solution

$${\tau }_x = \left({\tau }_2-{\tau }_1\right)\frac{\ln \frac{r_x}{r_1}}{\ln \frac{r_2}{r_1}}+{\tau }_1$$

(50)

with r₁ and r₂ being the internal and external radius of the tube, respectively; r_x the radius of an intermediate point between r₁ and r₂; and τ the temperature at the point considered.

The profile of the temperatures of a parametric section made in the FEMM model has been compared with the values obtained from complete analytical Equation (50), as can be seen in Figure 3. In this parametric section, the temperatures were measured for two concrete meshes, one with 2517 nodes and the other with 1,013,144 nodes.

The boundary conditions were the following. The highest temperature was on the outer surface of the cylinder (80 °C), and the lowest was on the inner face of the cylinder (30 °C). The top and bottom cylinder covers were considered perfect insulators. There was no heat source inside the cylinder. The initial and final conditions were the same, because it is a stationary process.

Figure 4 shows the results obtained for the two meshes analyzed with FEMM and the analytical solution for the parametric cut shown in Figure 3. We observed that the three curves were coincident, which allowed us to conclude that the solution with 2517 nodes was already a very good solution.

Table 1. Temperatures for the validation of FEMM.

Mesh [Nodes]	2517	1013144
Temperature (°C) a 0.75 m (Analytical ref: 72.7422 °C)	72.7585	72.7585
Difference of temperatures (°C)	0.0162	0.0162
Percentage difference (%)	0.0223	0.0223

shows the temperature values at a characteristic point located 0.75 m from the cylinder axis, and compares the temperature obtained by the exact analytical equation with the result obtained by FEMM for 2517 and 1,013,144 nodes. As it can be seen, the error converged for 2517, and was almost coincident with the analytical value, with a 0.0223% difference percentage.

Table 2. Metrics for the validation of FEMM vs. analytical solution. from the contrast of both results (the FEMM and studied model), we obtained the proposed metrics. R²—coefficient of determination; RMSPE—root mean square percentage error; MAEP—mean absolute percentage error; PBIAS—percentage bias.

Mesh [Nodes]	2517	1013144	References
R2 [0, +1]. Optimum: +1	1	1	[26]
RMSPE [−1, +1]. Optimum: 0	0.0001	0.0001	[27]
MAEP [−1, +1]. Optimum: 0	0.0001	0.0001	[27]
PBIAS [−1, +1]. Optimum: 0	0.0001	0.0001	[28]

shows the metrics used to validate the proposed mathematical models. Metrics are statistics used to assess the errors between the standard and proposed model. As stated above, our standard model is the FEMM. We compared the results obtained with the proposed models, and, from the contrast of both results (the FEMM and studied model), we obtained the proposed metrics. R² is the coefficient of determination. RMSPE is the root mean square percentage error. MAEP is the mean absolute percentage error. PBIAS is the percentage bias. The metrics indicate the validity of FEMM as a reference standard, because the values obtained are very close to the optimum.

The temperatures obtained were close enough to the analytical solution to consider the FEMM a good benchmark. In the distribution of errors (Figure 5), there is a bias respect to the zero error. This indicates a certain systemic error in the FEMM, which we assumed (0%, 0.03%; see Figure 5).

The simulations carried out were developed on a PC type Dell, Intel® Core (TM) i7-3820, 3.6 GHz, 32 GB RAM. The software, freely used, was Gmsh [29], as a CAD tool, mesher, and post-processing in 3D; FEMM as a CAD tool, meshing, processing, and post-processing using 2D FEM; and Dev C ++ [30] to develop the calculations in CM.

All of the numerical simulations that were carried out followed the same methodology [31], which is the following: formulation of the problem; geometric modelling and definition of the thermal domain; establishment of the boundary conditions and the initial conditions; generation of the mesh; simulation with processing of the results; and, finally, analysis of the results.

3.2. Numerical Simulation with $\left[M_{\lambda }\right]$

The model proposed for the validation consisted of a tube that is similar to the stator or the rotor of an electric machine, especially that of an asynchronous electric machine. It consisted of a tube of 2.00 m in height and 2.00 m in diameter. The wall has a thickness of 0.50 m. This wall is composed of two concentric tubes, with a thickness of 0.25 m for each one. The thermal conductivities of the inner and outer tubes were λ₁ = 50.0 Wm⁻¹K⁻¹ and λ₂ = 193.0 Wm⁻¹K⁻¹, respectively.

The wall domain, Ω_w, is defined as the volume between the outer surface and the inner surface of the tube. The enveloping domain, Ω_e, is defined as the external volume that surrounds the tube. The boundary conditions were the following: Ω_w is a thermal conductor, Ω_e is a thermal insulator, and there are no internal heat sources.

The boundary conditions were the following: the external temperature, T_{0_ext} = 30 °C and the internal temperature, T_{0_int} = 80 °C. It was considered that the heat advanced in a radial direction, that the upper and lower covers of the tube were perfect insulators and that there were no heat sources inside the tube. In addition, it was considered a pure heat transmission without convection or radiation. Therefore, it is a three-dimensional problem with axial symmetry.

To measure the temperatures, both in the 2D and 3D distributions, a parametric cut was made on a segment, in a radial direction, and there, the temperatures were measured, as can be seen in Figure 6.

The inner tube had a thermal conductivity of 50.0 Wm⁻¹K⁻¹, typical of steels, with a value between 46.6 and 51.9 Wm⁻¹K⁻¹. The outer tube had a thermal conductivity of 193.0 Wm⁻¹K⁻¹, typical of aluminum alloys, with a value between 95.3 and 222.0 Wm⁻¹K⁻¹.

In Figure 7, we analyzed a two-dimensional temperature distribution, which corresponds to the generatrix rectangle of revolution of the tube wall. This temperature distribution has been generated with FEMM using FEM. The temperature that was used as a reference for the validation of the results was that obtained by a horizontal parametric cut at half the height of the bi-dimensional distribution.

Then, in this section, we calculated, through $\left[M_{\lambda }\right]$, the heat conduction in the tube of Figure 6. The matrix $\left[M_{\lambda }\right]$ was proposed as an alternative to the matrices used in [18,19]. These results will be compared with those obtained with the new matrix $\left[M_{\tau }\right]$ in the Section 3.2.

On the one hand, a 2D FEMM calculation of 11,042 nodes was carried out, and on the other hand, five calculations were made with CM for an increasing number of nodes, ranging from 108 to 13,136 nodes, using $\left[M_{\lambda }\right]$. From 4461 nodes, the results were practically coincident between CM 3D and FEMM 2D, which we considered to be the solution to the problem, because being axisymmetric, its 3D equivalent mesh had many more nodes. The results obtained are shown in Figure 8.

The distribution of errors when comparing CM with the temperatures obtained in FEM is shown in Figure 9. There is a bias between the actual distribution of errors (bar chart in blue) and the theoretical distribution (Gauss bell type, red color chart). This bias is approximately 0.1%.

The errors, as expected, decreased with the increase in mesh density, as shown in

Table 3. Temperature at a point of the tube with two conductivities using $\left[M_{\lambda }\right]$.

Mesh [Nodes]	108	182	4461	9679	13,136
Temperature (°C) at 0.75 m (FEMM ref: 64.0653 °C)	67.8991	61.7719	64.0343	63.6987	63.6975
Difference of temperatures (°C)	3.8338	−2.2934	−0.0310	−0.3666	−0.3678
Percentage difference (%)	5.9842	−3.5798	−0.0484	−0.5722	−0.5741

.

The temperatures in Table 3 were taken at one point placed radially at 0.75 m from the center of the tube, and axially at half its height (1.00 m; see Figure 6).

Applying the following metrics, the determination coefficient (R²), the root mean square perceptual error (RMSPE), the mean absolute percentage error (MAEP), and percentage bias (PBIAS), the values obtained can be seen in

Table 4. Metrics for thermal simulation using $\left[M_{\lambda }\right]$.

Mesh [nodes]	108	182	4461	9679	13,136	References
R2 [0, +1]. Optimum: +1	0.9929	0.9982	1.0000	1.0000	1.0000	[26]
RMSPE [−1, +1]. Optimum: 0	0.0518	0.0122	0.0029	0.0013	0.0012	[27]
MAEP [−1, +1]. Optimum: 0	0.0471	0.0086	0.0025	0.0011	0.0011	[27]
PBIAS [−1, +1]. Optimum: 0	0.0450	−0.0051	0.0025	0.0011	0.0009	[28]

.

The R² values indicate a good adjustment of data in all of the comparatives. The RMSPE, MAEP, and PBIAS are beside the optimum. Even so, all of the indicators are in the optimal range. Looking at Table 2, we can assure that the biggest error committed is the RMSPE, whose value is 0.0518%. This error decreases as the node density increases. It is a more than acceptable value.

We compared a 3D CM model, with low density and tetrahedral meshes, with an axial symmetry model in 2D, solved by FEM, with a dense triangular mesh. The errors were not significant. The calculation process was greatly simplified. Remember that the number of the nodes indicates the number of equations and unknowns that solve the global matrix. This allows for us to see that there is no significant improvement when increasing the number of nodes, with denser meshes, as the results obtained with meshes of 9679 and 13,136 nodes were very similar, as far as errors are concerned.

3.3. Numerical Simulation with $\left[M_{\tau }\right]$

In the geometric model, the boundary conditions were the same as in the previous case. Therefore, the numerical simulation has been carried out with the tube indicated in Figure 6, and the FEMM model used is shown in Figure 7. Similarly, the three-dimensional parametric section can be seen in Figure 6. However, in this case, the thermal conductivity matrix was $\left[M_{\tau }\right]$.

The objective was to check the validity of the new thermal conductivity matrix $\left[M_{\tau }\right]$, and that it obtains better results than $\left[M_{\lambda }\right]$.

On the one hand, a calculation with a 2D FEMM of 11,042 nodes has been made, and, on the other hand, five calculations was made with CM for an increasing number of nodes ranging from 102 to 13,136 nodes, using $\left[M_{\tau }\right]$. From 6428 nodes, the results are practically coincident between 3D CM and 2D FEMM, which we consider as a solution to the problem. The results obtained are shown in the graph Figure 10.

The distribution of errors when comparing CM with the temperatures obtained in FEM is shown in Figure 11. As can be seen, the bias has decreased to 0.05%. Hence, operating with a new thermal conduction constitutive matrix $\left[M_{\tau }\right]$ produces a smaller error than operating with $\left[M_{\lambda }\right]$.

As we expected, the errors decrease as the mesh density increases, as can be seen in

Table 5. Temperature at a point of the tube with two conductivities using $\left[M_{\tau }\right]$.

Mesh [nodes]	102	108	2328	6428	12,895
1	68.4981	67.7963	63.7282	63.7047	63.6975
Difference of temperatures [°C]	4.8022	4.1004	0.0323	0.0088	0.0016
Percentage difference [%]	7.5393	6.4375	0.0507	0.0138	0.0025

.

The temperatures in Table 3 have been taken at one point placed radially at 0.75 m from the center of the tube, and axially at half its height (1.00 m).

Applying the metrics R², RMSPE, MAEP, and PBIAS, the following values were obtained (see

Table 6. Metrics for thermal simulation using $\left[M_{\tau }\right]$.

Mesh [nodes]	102	108	2328	6428	12,895	References
R2 [0, +1]. Optimum: +1	0.9910	0.9932	1.0000	1.0000	1.0000	[26]
RMSPE [−1, +1]. Optimum: 0	0.0576	0.0508	0.0026	0.0027	0.0008	[27]
MAEP [−1, +1]. Optimum: 0	0.0522	0.0462	0.0023	0.0024	0.0007	[27]
PBIAS [−1, +1]. Optimum: 0	0.0496	0.0442	0.0014	0.0022	0.0004	[28]

).

Observing Table 6, we conclude that the errors are lower when the mesh density is increased, and also that errors are very small for meshes with low densities. In any case, the errors quickly converge towards the optimum of the metrics.

Observing

Table 7. Temperature errors using $\left[M_{\lambda }\right]$ and $\left[M_{\tau }\right]$.

	FEMM ref: °C	Obtained temp. ºC	Error %
$\left[M_{\lambda }\right]$ 13,136 nodes	64.0653	63.6975	−0.5741
$\left[M_{\tau }\right]$ 12,895 nodes	63.6959	63.6975	0.0025

, it is verified that the thermal constitutive matrix $\left[M_{\tau }\right]$ that we provide in this article behaves much better than with $\left[M_{\lambda }\right]$.

Similarly, comparing the distribution of the error when we are using $\left[M_{\lambda }\right]$ (Figure 9) with the distribution of the error when we use $\left[M_{\tau }\right]$ (Figure 11), it can be seen that in the case of $\left[M_{\tau }\right]$, there is less bias, and it is closer to zero error when we are using $\left[M_{\lambda }\right]$.

The proposal of the new thermal constitutive matrix $\left[M_{\tau }\right]$ fundamentally differs in not using averaged values, but exact values such as the Cartesian coordinates of the tetrahedron nodes and the vectors of the dual surfacesm as it can be seen in Equations (47) and (48). This implies greater precision, as its calculation base is strictly geometric and not approximate, as it can be seen in the error histograms in Figure 9 and Figure 11, or if the values shown in Table 3, Table 4, Table 5 and Table 6 are compared. The matrix $\left[M_{\tau }\right]$ converges more evenly towards the solution than the matrix $\left[M_{\lambda }\right]$ does.

4. Conclusions

In this paper, a new constitutive matrix $\left[M_{\tau }\right]$ for thermal conduction for tetrahedral meshes, in a steady state thermal regime, has been developed through a new algebraic methodology, using the Cell Method. The results have been compared with those obtained for the same problem by means of the constitutive matrix $\left[M_{\lambda }\right]$, developed previously in other works. Taking a 2D axisymmetric model as reference, calculated with the finite element method, the errors obtained with the new matrix $\left[M_{\tau }\right]$ are of the order of 0.0025%, much lower than those obtained with $\left[M_{\lambda }\right]$. On the other hand, the finite formulation and its associated numerical method, the Cell Method, allows for great flexibility when mathematically modelling a physical phenomenon such as conduction heat transfer.

The main advantage of the method proposed in this article is its simplicity. The constitutive matrices developed by previous methods presented complex calculations, while the new constitutive matrix, proposed in this work, depends exclusively on the coordinates of the vertices of the tetrahedra, which constitutes the mesh.

In addition, the errors are much smaller with the new matrix, and this allows for meshes of smaller number of elements, obtaining the same precision with a lower temporal cost.

As it has been underlined before, the simplicity of the method and its greater precision mean that the new methodology can be applied to more complex problems, such as the calculation of the thermal heating of the rotor and the stator of an electrical machine of much more complex geometry, and with more complex physical properties and boundary conditions, even convective types. This is critical in the problems in transitory regime.