Computational Thermoelectricity Applied to Cooling Devices

Abstract

This chapter presents a numerical formulation within the finite element method in order to computationally simulate thermoelectric devices. For this purpose, a theoretical formulation based on nonequilibrium thermodynamics with historical notes is previously outlined. Then, a brief description of the finite element is reported to express the thermodynamics governing equations in an amenable form to be numerically discretized. Finally, several applications of cooling thermoelectrics are performed to highlight the benefits of the finite element method. In particular, a commercial thermoelectric device is simulated and several variables such as extracted heat, voltage drop, and temperature distributions inside the thermoelements are represented for different operating conditions. In conclusion, the present numerical tool could be used as a virtual laboratory for the design and optimization of Peltier cells.

1. Introduction

Thermoelectricity was discovered during the nineteenth century, but it has not been until the decade of 1950 when thermoelectric materials have experimented a quick expansion due to the techniques of doping [1]. Nowadays, the improvement of thermoelectric materials is based on the well-known figure-of-merit $Z$ , which depends on temperature, Seebeck coefficient and thermal and electric conductivities. A high value of $Z$ provides a good thermoelectric material for which the Seebeck coefficient and electrical conductivity are elevated and the thermal conductivity is low.

According to $Z$ , thermoelectric materials can be classified into two groups:

• Bulk materials such as Bi₂Te₃, PbTe, CoSb₃, SiGe, which present $Z\le 1$ [2]. In the last decades, a new range of bulk materials with $Z\le 1.5$ called nanostructured (for instance, phono-glass electron-crystal material) has emerged. Bulk materials are mainly used in high temperature applications: 700–900 K, [3].

• Low-dimensional materials, such as Sb₂Te₃ superlattices, are based on the quantum confinement effect of the electron charge carriers, which allows to obtain higher values of Seebeck coefficients. The main advantage of these materials is their figure-of-merit ( $Z\approx 1.5$ ); however, they are expensive and not adequate for mass production.

To sum up, despite the fact that values of $Z\approx 3$ have been reached for research applications [4], the best commercial thermoelectric material holds $Z\approx 1$ . Therefore, the material science community has a challenge to improve the $Z$ of thermoelectric materials and, consequently, to exploit the potential of these materials.

Regarding thermoelectric applications, there are two operating thermoelectric modes: coolers and generators, see Figure 1. Typically, all thermoelectric are composed of thermocouples that are made up of two semiconductors p- and n-type called thermoelements. Furthermore, tin weldings, copper bars for conducting the electric intensity $I$ and alumina for isolating the device are also present, as observed in Figure 1. The position of the hot end depends on the operation mode:

• Cooler (Figure 1 left), heat is extracted from the cold end by applying an electric current.

• Generator (Figure 1 right), an electric current is produced by the gradient of temperatures between cold and hot ends.

Figure 1.

Sketches of thermocouples operating as cooler (left) and generator (right).

Concerning cooling mode, the devices made out of thermoelectric materials have been restricted, among others, to four specific applications:

1. Refrigeration: Thermoelectric coolers present high reliability, no moving parts and reduction in size and weight. Nevertheless, the main drawback of thermoelectric refrigeration is that its efficiency is normally between 0.3 and 0.5 for a gradient of temperature of approximately 20°C, [5].

2. Electronic devices: These equipments, such as PC processors, generate considerable quantities of heat that cause troubles if the device junction temperature (approximately 85°C) is exceeded, [2]. Since conventional cooling techniques cannot completely fulfill the requirements of heat dissipation and due to the lack of space in electronic devices, the most adequate refrigeration option consists of using thermoelectric coolers combined with air/liquid cooling systems.

3. Automotive industry: Conventional techniques that use refrigerants such as R-134a contribute to global warming. In this connection, thermoelectric coolers could solve the problem. Furthermore, studies show that acceptable efficiency in the range of 0.4–0.8 can be achieved under ambient temperature from 46 to 30°C in the case of a truck cab, [6].

4. Air-conditioning: In this field, thermoelectric air-conditioning is not still as efficient as the vapor-compression system. For instance, efficiency of vapor-compression and thermoelectric air-conditioning systems are in the range of 2.6–3.0 and 0.38–0.45, respectively, see [7].

With regard to thermoelectric generator, there exists a great interest in the use of generator devices made out of thermoelectric materials to produce energy from residual sources such as waste heat. These applications are commonly denominated energy harvesting and, nowadays, have a big potential. For this reason, most of the research works focus on thermoelectric generators.

The aim of the present chapter is to develop a numerical formulation within the finite element (FE) method to computationally model thermoelectric devices. Previously, a review of the different numerical techniques that are applied to model thermoelectricity is reported. Then, the FE formulation is introduced and several applications are developed in order to highlight the main advantages of the application of the FE method to thermoelectricity. Notice that this chapter is only focused on cooling applications for the sake of brevity.

This chapter is structured as follows: Firstly, Section 2 presents a state of art on the different numerical techniques to simulate thermoelectricity. Secondly, a thermodynamically consistent formulation within the nonequilibrium thermodynamic formalism is developed in Section 3 to obtain the governing equations, which are composed of balance laws, transport equations and boundary conditions. Thirdly, Section 4 applies the FE method to rewrite the governing equations in an amenable form to be computationally solved. Finally, several examples such as convergence tests, comparisons with analytical solutions and miniaturizations are conducted in Section 5 and the conclusions are summarized in Section 6.

For the sake of clarity, the notation used in the present chapter is shown in Table 1.

Symbol	Description
	Theoretical formulation
$\mathrm{\Omega }$ , ${\mathrm{\Omega }}^{\infty }$ , $\mathrm{\Gamma }$ , $\mathbf{n}$	Domain, surrounding, boundary and outward normal
$\rho =\left\{{\rho }_q^f,{\rho }_m\right\}$	Free electric charge and mass density
$\mathbf{J}=\left\{\mathit{j},{\mathit{j}}_s\right\}$	Electric and entropy fluxes
$\mathrm{\Xi }={\mathrm{\Sigma }}_{\mathrm{s}}$	Entropy production
$\mathrm{d}$ , $\delta$	Exact differential and variation
$t$ , $\mathrm{s}$ , $T$ , $V$	time, entropy, temperature and voltage
$\mathit{q},\mathit{x}$ , $$	Heat flux and Cartesian coordinates
$c$ , $\alpha$ , $\kappa$ , $\gamma$	Heat capacity, Seebeck coefficient, and thermal and electric conductivities
	Finite element formulation
$g_a=\left\{{\tilde{V}}_a,{\tilde{T}}_a,{\dot{\tilde{T}}}_a\right\}$	Degrees of freedom
${\mathrm{\Omega }}_e$ , ${\mathrm{\Gamma }}_e$	Subdomain and boundary of any finite element
$\mathcal{R}$ , $\mathcal{K}$ , $\mathcal{C}$	Residual, tangent stiffness and capacity matrices
$a$ , $b$	Global numbering of elements
$k$	Newton-Raphson counter
$c_1$ , $c_2$	Time integration parameters
	Applications
$T_c$ , $T_h$	Cold and hot temperatures
$Q_c$ , $Q_h$	Extracted heat at cold and hot ends
$I$ , $Z$	Intensity and figure-of-merit

Table 1.

Table of symbols.

2. Review of numerical methods applied to thermoelectricity

Despite the fact that there exists a large amount of works that use simplified analytical models [8, 9, 10] and semi-analytical techniques such as Laplace transforms [11], numerical modeling of thermoelectricity is a fundamental aspect in order to design and optimize sophisticated devices made out of thermoelectric materials. In this connection, there are at least five numerical methods applied to solve the set of two coupled partial differential equations (PDE) of thermoelectricity. For more details, readers are referred to [12], in which some of these numerical methods are compared and discussed.

• Network or circuit simulator method, which can be considered as a representation of the natural convection and uses the Kirchhoff’s theorems. The main characteristic of this method is the discretization of the spatial coordinates while time is assumed to be a continuous variable. For instance, the method developed in [13] considers the Thomson effect, which is often neglected in the models for simplicity, and also the temperature-dependency of all material parameters.

• Finite difference method is the simplest numerical technique and approximates the set of PDE by using the Taylor’s theorem. An example of application can be found in [14]. In this work, the authors study the production of electricity by using the waste heat of a combustion chamber. The temperature and mass flow of the flue gases and the load resistance, which are key parameters for the generation, are included in the model, as well as the elements of the generators (heat exchangers, ceramics, unions) and all the thermoelectric phenomena. As a main novelty of the model, the temperature drop of the flue gases while circulating along the thermoelectric generator is also included. As shown by the authors, this thermoelectric model based on finite differences predicts the thermoelectric generation with a relative error within ±12%. Also in [15] a computational model for thermoelectric self-cooling, based on the implicit finite differences, is presented. The model is capable of simulating both the steady and the transient state of the system and considers all the thermoelectric effects (Seebeck, Thomson, Peltier and Joule) including temperature-dependency of all properties, and exhibits a good correspondence with the experimental measurements (±12% of deviation).

• Boundary element method [16] solves the PDE by integral equations and is applied to problems for which Green’s functions can be calculated. The main shortcoming of this method is its difficulty for solving nonlinear problems and, for this reason, there are few works in the literature applied to thermoelectricity. For instance, this method is used in [17] to study the two-dimensional temperature distributions in thermoelectric devices.

• Finite volume method is also used to evaluate the set of coupled PDE of thermoelectricity in [18]. In this case, the authors numerically study the thermoelectric performance of composite thermoelectric devices, in terms of produced electrical current, Ohmic and Seebeck potentials, power output, heat input and conversion efficiency for various hot surface temperatures, load resistance, semiconductor thickness and convection heat transfer coefficient.

• FE method [19] is the most advanced technique for solving PDE in science and engineering. Furthermore, this method is optimal for applying to nonlinear coupled problems as thermoelectricity. In this connection, the authors of the present chapter have extensive experience in developing thermodynamic consistent FE formulations, which are implemented in the research code FEAP [20], see [21, 22, 23, 24, 25] and [26]. Regarding commercial FE codes, the authors in [27, 28] use ANSYS to perform a comprehensive numerical analysis focusing on the cooling performance of miniaturized thermoelectric coolers for microelectronics applications. In particular, they study the effects of parameters such as the load current, geometric size. Finally, the software COMSOL is used in [29] to study the thermal stresses in Peltier cells.

In short, the main disadvantage of the boundary element is, as commented, its intrinsic difficulty to solve nonlinear problems. On the contrary, all the other methods can be applied to solve nonlinear thermoelectricity. The FE and finite volume methods require of complex mathematical developments and they are commonly used by the continuum mechanics community. In contrast, the finite difference method is the most direct approach to numerically solve PDE. Finally, the circuit method is preferably used by the electric engineering community.

Despite the difficulty of the FE method, it is used in the present chapter since it presents two main advantages: i) It is a very general an efficient method, and ii) it is easy to improve the accuracy of the solutions by two approaches: mesh refinement and/or increasing the order of the interpolation functions.

3. Governing equations

Thermoelectricity couples electric and thermal energies by five separated effects: Fourier, Ohm, Seebeck, Thomson and Joule (see historical notes 1 and 2). All of them are thermal and electric fluxes that must satisfy the balances of electric charge and of entropy. The aim of this section is to state both balances and, in addition, to obtain the transport equations that couple electric and thermal energies. Finally, this section reports the boundary conditions that are required for a proper formulation of the set of coupled PDE that governs thermoelectricity.

Historical note 1

• Jean-Baptiste Joseph Fourier (1768–1830) was a French mathematician and physicist. Coming of a humble family, he showed such proficiency in mathematics in his early years that he began his career in 1790 as a teacher of mathematics in the school he first attended, the École Royale Militaire of Auxerre. Fourier also joined the local revolutionary committee in his own district promoting the ideals of the French Revolution, and his life was in danger several times for this reason. However, he is best remembered for his work on flow of heat, through his Théorie analytique de la chaleur (1872). This work is the basis of what is today known as the Fourier transform.

• Georg Ohm (1789–1854) was a German physicist who summarized the most relevant aspects of Ohm’s law in his work Die galvanische kette, mathematish bearbeitet (1827). At that time he was professor of mathematics at the Jesuits College at Cologne, where he had a well-equipped laboratory in which he started his first experiments with electricity. Despite the great influence of his work in the theory and applications of electric currents, it was not before 1841 when he started to be recognized, being awarded by the Royal Society of London. Also in his honor, the physical unit measuring the electrical resistance ohm, was named for him after his death.

• James Prescott Joule (1818–1889) was an English physicist best known by his paper On the production of heat by voltaic electricity (1840). In this work, he stated the principles of the law that has his name, Joules law: mathematical description of the rate at which the resistance of a circuit converts electric energy into heat. It was also Joule who formed the basis of the first law of thermodynamics, since he established that the various forms of energy (mechanical, electrical and heat) were in essence the same and could be changed. Other relevant contribution emerged in 1852 when Joule and Thomson discovered the Joule-Thomson effect, which is commonly exploited in thermal machines and played a crucial role in the development of the refrigeration industry in the nineteenth century.

3.1. Balance of electric charges and entropy

The term continuum physics refers to several branches of physics such electromagnetism and nonequilibrium thermodynamic for which the matter consists of material points that are localized by the Euclidean position vector $\mathit{x}$ . This fact is mathematically grounded on the continuum hypothesis [30], which allows to work with balance equations that state the conservation, production or annihilation of certain quantities such electric charge and specific entropy.

Consider a closed system of domain $\mathrm{\Omega }$ , boundary $\mathrm{\Gamma }$ with outward normal $\mathbf{n}$ and its surrounding ${\mathrm{\Omega }}^{\infty }$ . A general expression of a balance equation for the quantity $\mathrm{\Upsilon }$ at time $t$ and in any material point $\mathit{x}$ is given by:

$$\frac{d}{dt}{\int }_{\mathrm{\Omega }}\rho \left(\mathrm{\Upsilon };\mathit{x},t\right)\mathrm{d\Omega }=-{\int }_{\mathrm{\Omega }}\nabla \cdot \mathbf{J}\left(\mathrm{\Upsilon };\mathit{x},t\right)\mathrm{d\Omega }+{\int }_{\mathrm{\Omega }}\mathrm{\Xi }\left(\mathrm{\Upsilon };\mathit{x},t\right)\mathrm{d\Omega },$$ E1

where $\rho$ , $\mathbf{J}$ and $\mathrm{\Xi }$ denote the $\mathrm{\Upsilon }$ -density, $\mathrm{\Upsilon }$ -current density and $\mathrm{\Upsilon }$ -production, respectively. This equation balances the quantity $\mathrm{\Upsilon }$ since it contains information on:

• The quantity of $\mathrm{\Upsilon }$ inside $\mathrm{\Omega }$ , term on the left-hand side.

• Inflow quantity from ${\mathrm{\Omega }}^{\infty }$ to $\mathrm{\Omega }$ through the boundary $\mathrm{\Gamma }$ , first term on the right-hand side (after application of divergence theorem).

• Production/annihilation of $\mathrm{\Upsilon }$ inside $\mathrm{\Omega }$ , second term on the right-hand side.

In thermoelectricity, both electric and thermal fields are present and, consequently, there are two $\mathrm{\Upsilon }$ quantities as shown in Table 2. As observed, the quantities are the free electric charge ${\rho }_q^f$ and the specific entropy $\mathrm{s}$ ; the current densities are the electric $\mathit{j}$ and entropy ${\mathit{j}}_{\mathrm{s}}$ fluxes. Regarding the production terms, there is not production for the electric field due to the fact that the electric charge can neither be created nor destroyed. Consequently, the electric balance becomes a conservation equation. On the contrary, the entropy production ${\mathrm{\Sigma }}_{\mathrm{s}}$ is due to the irreversibilities (also called dissipations) generated by both the Fourier heat transport and the Joule heating.

Quantity	$\mathrm{\Upsilon }$	Current density	$\mathbf{J}$	Production	$\mathrm{\Xi }$
Free electric charge	${\rho }_q^f$	Electric flux	$\mathit{j}$	–	–
Specific entropy	${\rho }_m\mathrm{s}$	Entropy flux	${\mathit{j}}_{\mathrm{s}}$	Entropy production	${\mathrm{\Sigma }}_{\mathrm{s}}$

Table 2.

Quantities to be balanced in thermoelectricity.

In the present chapter, it is assumed that there are not free electric charges; this is a first and good approximation for thermoelectric devices made out of semiconductors. Under this assumption and introducing the quantities of Table 2 in (1), the thermoelectric balance equations in local forms read:

$$0=\nabla \cdot \mathit{j},{\rho }_m\dot{\mathrm{s}}=-\nabla \cdot {\mathit{j}}_{\mathrm{s}}-\frac{\mathit{j}\cdot \nabla V}{T},$$ E2

where the second term on the right-hand side of the right equation represents the Joule heating and $V$ , $T$ denote voltage and temperature, respectively.

The entropy balance of (2) can be manipulated by considering that ${\mathit{j}}_s=\mathit{q}/T$ , where $\mathit{q}$ is the heat flux to give:

$${\rho }_m\dot{\mathrm{s}}=-\nabla \cdot \left(\frac{\mathit{q}}{T}\right)\underset{{\mathrm{\Sigma }}_{\mathrm{s}}}{\underbrace{-\frac{\mathit{q}\cdot \nabla T}{T^2}-\frac{\mathit{j}\cdot \nabla V}{T}}}.$$ E3

Now, the classical heat balance equations can be obtained by manipulating (3) and by taking into account that $\dot{\mathrm{s}}=c\dot{T}/T$ , where $c$ is the specific heat capacity, to give:

$${\rho }_{m}c\dot{T}=-\nabla \cdot \mathit{q}-\mathit{j}\cdot \nabla V.$$ E4

As observed, the Joule heating is a heat source since the electric energy is irreversibly converted into heat. In addition, this term is nonlinear due to its quadratic dependency of the $\nabla V$ term, as can be observed from (6).

3.2. Transport equations

Within the nonequilibrium thermodynamic formalism [31], the transport equations are obtained from the entropy production ${\mathrm{\Sigma }}_{\mathrm{s}}$ of (3). Concretely, in a first and good approximation the entropy production can be expressed as the sum of products of fluxes and driving forces. In matrix form:

$$\left\{\begin{array}{l}-\frac{\nabla V}{T}\\ \\ -\frac{\nabla T}{T^2}\end{array}\right\}=\left[\begin{array}{ll}\frac{{\alpha }^2}{\kappa }+\frac{1}{\mathit{\gamma T}}& -\frac{\alpha }{\mathit{\kappa T}}\\ & \\ -\frac{\alpha }{\mathit{\kappa T}}& \frac{1}{\kappa T^2}\end{array}\right]\left\{\begin{array}{l}\\ \begin{array}{l}\mathit{j}\\ \\ \mathit{q}\end{array}\\ \end{array}\right\},$$ E5

where $\gamma$ , $\kappa$ and $\alpha$ denote electric and thermal conductivities and Seebeck coefficient, respectively. Finally, the set of two coupled equations of thermoelectricity becomes:

$$\mathit{j}=-\gamma \left(T\right)\nabla V-\alpha \left(T\right)\gamma \left(T\right)\nabla T,\mathit{q}=-\kappa \left(T\right)\nabla T+\alpha \left(T\right)T\mathit{j}.$$ E6

Historical note 2

The thermoelectric coupling was discovered by Thomas Johann Seebeck (1770–1831), Athanase Peltier (1785–1845) and William Thomson (1824–1907). The former discovered the Seebeck effect in a casual experiment; he initially believed that the deflection experienced by a compass magnet in a circuit made from two different metals with junctions exposed to a temperature gradient was due to magnetism, instead of considering that an electric current was induced by the temperature difference. For that reason, for a century it had no application. The French physicist Peltier realized that the passage of an electric current would induce heating or cooling at the junction of two dissimilar metals. Today, this effect is known called Peltier-Seebeck effect and is the physical basis of thermocouples. Finally, Thomson, later Lord Kelvin, showed the appearance of a temperature gradient when current flows in a material, resulting in a cooling or heating effect of a different nature from Joule effect. Thomson is also best known for giving a comprehensive explanation of the Seebeck and Peltier effects and describing their interrelationships, known as Kelvin relations. The definition of the absolute temperature scale, based on the degree Celsius, was also other of his achievements.

In semiconductors, the material properties $\gamma$ , $\kappa$ and $\alpha$ commonly depend on temperature. For instance, for the metal-metalloid alloy ${\mathrm{Bi}}_2{\mathrm{Te}}_3$ the dependency on temperature of these properties is shown in Figure 2. As observed, both $\gamma$ and $\kappa$ exhibit a strong temperature dependency that results in a material nonlinearity. Explicitly, the temperature dependency of these material properties can be fitted to second order polynomials to give:

$$\begin{array}{llll}\gamma \left(T\right)={\gamma }_0+{\gamma }_{1}T+{\gamma }_2{T}^2,& \kappa \left(T\right)={\kappa }_0+{\kappa }_{1}T+{\kappa }_2{T}^2,& \alpha \left(T\right)={\alpha }_0+{\alpha }_{1}T+{\alpha }_2{T}^2.& \end{array}$$ E7

Figure 2.

Temperature-dependency of material properties for the ${\mathrm{Bi}}_2{\mathrm{Te}}_3$ material. Picture taken from [21].

Notice that the temperature-dependency of $\alpha \left(T\right)$ is the responsible of the Thomson’s effect and, consequently, the inclusion of material nonlinearities is a requirement for a proper modeling of thermoelectricity.

3.3. Boundary conditions

The boundary conditions are additional constrains at the boundary $\mathrm{\Gamma }$ that are present at all the boundary value problems. These conditions guarantee the existence of a unique solution and, therefore, the problem is well posed. There exist two types of boundary conditions:

• Dirichlet conditions (also called first-type or essential) specify the value of the degrees of freedom at the boundary.

• Neumann conditions (second-type or natural) give the value of the normal derivative of the fluxes at the boundary.

In thermoelectricity, these conditions read:

$$\begin{array}{rrr}\text{Dirichlet}\text{type}:& V=\overline{V},& T=\overline{T},\\ \text{Neumann}\text{type}:& \mathit{j}\cdot \mathbf{n}=\overline{j},& \mathit{q}\cdot \mathbf{n}=\overline{q},\end{array}$$ E8

where $\overline{V}$ , $\overline{T}$ , $\overline{j}$ and $\overline{q}$ denote prescribed voltage, temperature, electric flux and thermal flux, respectively.

4. Finite element formulation

As commented, the FE is probably the most advanced method for the solution of coupled PDE, but it involves complex mathematical concepts. For the sake of clarity and in order to present a comprehensible FE formulation, this section reports the FE formulation of thermoelectricity by using the following steps:

Historical note 3

The FE method was developed in the 1960s, among others, by the Greek John Argyris (1913–2004) at the University of Sttutgart, by the American Ray William Clough (1920–2016) at the University of California (Berkeley) and by the Anglo-Polish Olgierd Zienkiewicz (1921–2009) at the University of Swansea. The main contribution of the latter was to recognize the general potential of FE to resolve problems in areas outside solid mechanics.

i) The continuum domain $\mathrm{\Omega }$ is divided into subdomains called finite elements ${\mathrm{\Omega }}_e$ interconnected at nodal points $a$ , see Figure 3. Mathematically, this discretization is expressed by $\mathrm{\Omega }\approx {\cup }_{e=1}^{n_{\mathit{el}}}{\mathrm{\Omega }}_{e_i}$ . There exists a wide amount of finite element types depending on their number of nodes, spatial dimensions, etc. For more details, readers are referred to [19]. In the present work, three-dimensional eight noded elements, as that shown in Figure 3, are considered.

Figure 3.

Finite element discretization of a continuum domain and representation of an eight noded element.

ii) The values at nodes are called degrees of freedom $\mathit{g}$ and are the problem unknowns. In thermoelectricity, they are voltage $V$ , temperature $T$ and their time derivatives:

$${\mathit{g}}_a=\left\{{\tilde{V}}_a,{\tilde{T}}_a,{\dot{\tilde{T}}}_a\right\}.$$ E9

iii) A set of functions denominated “shape functions” $\mathcal{N}$ are chosen to interpolate the unknowns and the spatial dimensions within each finite element in terms of their nodal values. Obviously, these functions must satisfy several mathematical requirements that, for the sake of brevity, are not reported in the present work. Furthermore, there are many types of shape functions depending on their polynomial type, degree of interpolation, etc. In this chapter, standard shape functions of Lagragian type, which are linear functions of the degrees of freedom, are used and the numerical interpolations become:

$$\mathit{x}\approx {\mathcal{N}}_a{\tilde{\mathit{x}}}_a,V\approx {\mathcal{N}}_a{\tilde{V}}_a,T\approx {\mathcal{N}}_a{\tilde{T}}_a,\dot{T}\approx {\mathcal{N}}_a{\dot{\tilde{T}}}_a,$$ E10

where the Einstein summation convention is used.

iv) The balance and constitutive equations and the boundary conditions reported in Section 3 represent the strong form of the thermoelectric problem; namely, the coupled PDE depend on the second spatial derivatives of $V$ and $T$ . The FE uses “weakened” forms for which $V$ and $T$ are first spatial derivatives. For this purpose, the governing equations must be manipulated by using several approaches such as the principle of virtual work, the weighted residual, Hamilton’s principle, etc. In the present work, the weighted residual method is used since it is probably the most systematic approach; it consists of three steps:

a) The balance equations of (2) (left) and (4) are multiplied by test functions of the degrees of freedom $\mathit{\delta V}$ and $\mathit{\delta T}$ .

b) The divergence theorem is applied to both arising equation.

c) The Neumann boundary conditions are enforced and the weak forms yield:

$$\begin{array}{lll}{\int }_{\mathrm{\Omega }}\nabla \mathit{\delta V}\cdot \mathit{j}\mathrm{d\Omega }-{\oint }_{\mathrm{\Gamma }}\mathit{\delta V}\overline{j}\mathrm{d\Gamma }& =& 0,\\ {\int }_{\mathrm{\Omega }}\left[\nabla \mathit{\delta T}\cdot \mathit{q}+\mathit{\delta T}\left(\mathit{j}\cdot \nabla V+{\rho }_{m}c\dot{T}\right)\right]\mathrm{d\Omega }-{\oint }_{\mathrm{\Gamma }}\mathit{\delta T}\overline{q}\mathrm{d\Gamma }& =& 0.\end{array}$$ E11

v) In this step, similar discretization to those of (10) are applied to the test functions and their spatial derivatives:

$$\begin{array}{ll}\mathit{\delta V}\approx {\mathcal{N}}_a\delta {\tilde{V}}_a,& \mathit{\delta T}\approx {\mathcal{N}}_a\delta {\tilde{T}}_a,\\ \nabla \mathit{\delta V}\approx \nabla {\mathcal{N}}_a\delta {\tilde{V}}_a={\mathcal{B}}_a\delta {\tilde{V}}_a,& \nabla \mathit{\delta T}\approx \nabla {\mathcal{N}}_a\delta {\tilde{T}}_a={\mathcal{B}}_a\delta {\tilde{T}}_a.\end{array}$$ E12

Introducing these expressions in the weak form of (11), the FE residuals for each finite element of domain ${\mathrm{\Omega }}_e$ of boundary ${\mathrm{\Gamma }}_e$ become:

$$\begin{array}{lll}{\mathcal{R}}_b^V& =& {\int }_{{\mathrm{\Omega }}_e}{\mathcal{B}}_b^{\mathrm{{\rm T}}}\mathit{j}\mathrm{d}{\mathrm{\Omega }}_e-{\oint }_{{\mathrm{\Gamma }}_e}{\mathcal{N}}_b\overline{j}\mathrm{d}{\mathrm{\Gamma }}_e,\\ {\mathcal{R}}_b^T& =& {\int }_{{\mathrm{\Omega }}_e}\left[{\mathcal{B}}_b^{\mathrm{{\rm T}}}+{\mathcal{N}}_b\left(\mathit{j}\cdot \nabla V+{\rho }_{m}c\dot{T}\right)\right]\mathrm{d}{\mathrm{\Omega }}_e-{\oint }_{{\mathrm{\Gamma }}_e}{\mathcal{N}}_b\overline{q}\mathrm{d}{\mathrm{\Gamma }}_e,\end{array}$$ E13

vi) Finally, the solution is calculated by solving a set of nonlinear transient equations that involves three steps:

(a) The time interval is divided into small time increments $\mathrm{\Delta }t$ .

(b) The time derivatives are replaced by discrete forms using time integration techniques such as backward or forward finite differences, see [19].

(c) Since the thermoelectric problem is nonlinear due to both the temperature-dependency of the material properties and the Joule heating, the resulting nonlinear algebraic problem is linearized at each time increment by using a numerical algorithm. For instance, the Newton-Raphson algorithm use $k$ iterations for the linearization of the residuals:

$${\left.{\mathcal{R}}_a\right|}^k=-{\left.\frac{\partial {\mathcal{R}}_a}{\partial g_b}\right|}^k\mathrm{d}{g}_b^k\Rightarrow -{\left.\frac{\partial {\mathcal{R}}_a}{\partial g_b}\right|}^k=c_1{\mathcal{K}}_{\mathit{ab}}+c_2{\mathcal{C}}_{\mathit{ab}},$$ E14

where $a$ , $b$ refer the local numbering of two generic FE nodes, the parameters $c_1$ and $c_2$ depend on the time integration scheme, see [19], and the consistent tangent matrices $\mathcal{K}$ and the nonzero capacity matrix $\mathcal{C}$ are derived at each iteration as:

$$\begin{array}{l}{\mathcal{K}}_{\mathit{ab}}^{\mathit{VV}}=-\frac{\partial {\mathcal{R}}_a^V}{\partial {\tilde{V}}_b}=-{\int }_{{\mathrm{\Omega }}_e}{\mathcal{B}}_a^{\mathrm{{\rm T}}}\frac{\partial \mathit{j}}{\partial {\tilde{V}}_b}\mathrm{d}{\mathrm{\Omega }}_e,\\ {\mathcal{K}}_{\mathit{ab}}^{\mathit{VT}}=-\frac{\partial {\mathcal{R}}_a^V}{\partial {\tilde{T}}_b}=-{\int }_{{\mathrm{\Omega }}_e}{\mathcal{B}}_a^{\mathrm{{\rm T}}}\frac{\partial \mathit{j}}{\partial {\tilde{T}}_b}\mathrm{d}{\mathrm{\Omega }}_e,\\ {\mathcal{K}}_{\mathit{ab}}^{\mathit{TV}}=-\frac{\partial {\mathcal{R}}_a^T}{\partial {\tilde{V}}_b}=-{\int }_{{\mathrm{\Omega }}_e}\left[{\mathcal{B}}_a^{\mathrm{{\rm T}}}\frac{\partial \mathit{q}}{\partial {\tilde{V}}_b}+{\mathcal{N}}_a\frac{\partial \left(\mathit{j}\cdot \nabla V\right)}{\partial {\tilde{V}}_b}\right]\mathrm{d}{\mathrm{\Omega }}_e,\\ {\mathcal{K}}_{\mathit{ab}}^{\mathit{TT}}=-\frac{\partial {\mathcal{R}}_a^T}{\partial {\tilde{T}}_b}=-{\int }_{{\mathrm{\Omega }}_e}\left[{\mathcal{B}}_a^{{\rm T}}\frac{\partial \mathit{q}}{\partial {\tilde{T}}_b}+{\mathcal{N}}_a\frac{\partial \left(\mathit{j}\cdot \nabla V\right)}{\partial {\tilde{V}}_b}\right]\mathrm{d}{\mathrm{\Omega }}_e,\\ {\mathcal{C}}_{\mathit{ab}}^{\mathit{TT}}=-\frac{\partial {\mathcal{R}}_a^T}{\partial {\dot{\tilde{T}}}_b}=-{\int }_{{\mathrm{\Omega }}_e}{\mathcal{N}}_a{\rho }_{m}c{\mathcal{N}}_b\mathrm{d}{\mathrm{\Omega }}_e,\end{array}$$ E15

where the derivatives are reported in [24] and are not repeated here for the sake of brevity.

Finally, the solution is updated by using ${\mathit{g}}_b^{k+1}={\mathit{g}}_b^k+d{\mathit{g}}_b^k$ until the convergence is reached. In particular, the Newton-Raphson algorithm exhibits a quadratic asymptotic rate of convergence. The present FE formulation is coded in the research software FEAP [20], which belongs to the University of California at Berkeley (USA).

5. Numerical examples

As commented, the present numerical tool can be used as a “virtual laboratory” that allows to numerically design and optimize devices made out of thermoelectric materials. In this connection, the aim of this section is to show several features of the FE code. In particular, a commercial Peltier cell for cooling applications is simulated and several variables such as heat extracted, voltage and temperature distributions are analyzed.

A CP1.4-127-045 thermoelectric cooling module composed of 127 thermocouples electrically connected in series, as observed in Figure 4 (left), and manufactured by LairdTech [32] is modeled. The technical specifications of this device are: maximum intensity $I=8.7$ (A), maximum extracted heat $Q_c=82.01$ (W) with voltage drop $V$ = 15.33 (V) and under $T_h=T_c={50}^{\circ }$ C, where $T_c$ and $T_h$ refer to the temperature at cold and hot ends, respectively.

Figure 4.

Left: Sketch of the thermoelectric cooling module. Right: Finite element mesh and boundary conditions applied to half of the thermocouple.

Numerically, only half of the thermocouple needs to be modeled due to its symmetry and, consequently, the computational time is reduced; the FE mesh composed of 12,670 eight noded elements is shown in Figure 4 (right). On the one hand, the Neumann boundary conditions, namely, the electric current is enforced by the specific two-dimensional FE developed in [21]. This element can be also used to take into account convection and radiation phenomena among thermocouples; however, both phenomena are ignored in the present chapter. On the other hand, the Dirichlet boundary conditions are applied by setting the voltage reference $V=0$ and by prescribing temperatures at cold and hot ends, see Figure 4 (right). Finally, the coefficients of the temperature-dependent material properties of (7) are also reported in [21].

5.1. Convergence tests

In order to study the accuracy of the FE simulation and since the problem is nonlinear, convergence tests must be carried out. In particular, this section presents two tests: h and residual convergences. The former refers to the improvement of the solution with decreasing the size of the finite elements—this convergence is called h-refinement by the computational mechanics community—and the latter guarantees a proper performance of the Newton-Raphson algorithm to solve the nonlinearities.

For this purpose, the most detrimental operating condition is simulated: the boundary conditions are $T_c=T_h={50}^{\circ }$ C and $I=8.7$ (A). Figure 5 shows the h-refinement (left) and residual convergence (right) for the potential drop $V$ (top) and the extracted heat at the cold end $Q_c$ (bottom).

Figure 5.

Left: mesh convergence, voltage drop (top) and heat extracted (bottom) vs. number of elements. Right: residual convergence, residual norm of voltage (top) and heat extracted (bottom) vs. number of iterations.

As observed in the left figure, a mesh composed of 12,670 elements gives a perfect convergence for both electrical and thermal variables. The potential drop converges faster than the extracted heat; due to the fact that the voltage is directly calculated while the heat depends on spatial derivatives. Regarding the residual convergence and as observed in the right figure, both electrical and thermal variables require three iterations of the Newton-Raphson algorithm (the residual norm becomes ${10}^{-16}$ ) to minimize the residual due to the relative nonlinearity of the problem. Notice that each iteration requires 0.46 (s) in a 1.6 (GHz) Intel Core i5 processor.

5.2. Comparisons of extracted heat and voltage drop

The aim of this section is to compare analytical and FE solutions for the commercial cell under $I=1.7$ (A) and $T_h={50}^{\circ }$ C. The analytical solutions are obtained from [33] and are typically used by designers. For this aim, Figure 6 shows the voltage drop (left) and the extracted heat (right) versus the temperature at the cold end.

Figure 6.

Voltage drop (left) and extracted heat (right) vs. temperature at the cold end for $I=1.7$ (A), $T_h={50}^{\circ }$ C.

As observed, analytical and numerical solutions for the voltage disagree due to the electrical loss at the corners of the thermoelements. This drop loss is captured by the FE but not by the one-dimensional analytical solution, as was concluded in [21] by comparing analytical, numerical and experimental solutions. Therefore, for a proper design of thermoelectric cells the numerical tool is more appropriated than the simple analytical solutions. For the extracted heat, despite the fact that the slope is slightly different, both analytical and numerical solutions agree.

Regarding the drop loss at the corners of the thermoelement, Figure 7 shows contour plots of the electric currents along horizontal (left) and vertical (right) axes. As observed, the fluxes produce vortices-like phenomena and, consequently, part of the prescribed intensity does not go to the thermoelement.

Figure 7.

Contour plots of electric current along horizontal (left) and vertical (right) axes for $I=1.7$ (A), $T_h=T_c={50}^{\circ }$ C.

Another important aspect to be considered for the proper operation of Peltier cells is the maximum temperature inside the thermoelements since an overheating produces thermal stresses that could break the semiconductors from a mechanical point of view. For that matter, Figure 8 shows the maximum and minimum temperature inside the thermoelements for two operating conditions: $I=1.7$ and $I=3$ (A), at $T_h={50}^{\circ }$ C.

Figure 8.

Minimum and maximum temperature inside thermoelements vs. temperature at the cold end for $I=1.7$ (left) and $I=3$ (A) (right), $T_h={50}^{\circ }$ C.

As observed, for $I=3$ (A) the maximum temperature becomes approximately 80 ${}^{\circ }$ C due to the increasing of the Joule heating that depends on the prescribed intensity. To sum up, the present numerical tool can be used for a proper design and optimization of cooling devices from both thermoelectric and thermomechanic interactions.

5.3. Miniaturization of thermoelements

The last example deals with the minimization of thermoelements. Currently, there exists a trend towards minimization of thermoelectric devices for two reasons: decreasing power consumption and reducing the size of the devices.

Figure 9 shows the voltage drop (left) and the extracted heat (right) versus the thermoelement length for the detrimental case $I=8.7$ (A) and $T_c=T_h={50}^{\circ }$ C. As observed, decreasing the size of thermocouple implies to reduce the heat extracted and increase the potential drop. Obviously, both results are detrimental since the coefficient of performance of the cell becomes worse. This fact is due to the Joule heating that is a volumetric effect and, consequently, decreasing the size and maintaining constant the prescribed intensity results in overheating of the thermoelements. Again, the present modeling tool could be used for calibrating the miniaturization of devices by setting the applied intensity in order to achieve an efficient functioning.

Figure 9.

Voltage drop (left) and extracted heat (right) vs. the thermoelement length.

6. Concluding remarks

This chapter has presented a thermodynamically consistent formulation to obtain the governing equations of thermoelectricity, from a theoretical point of view. Then, a nonlinear numerical formulation within the finite element method is developed. The nonlinearities emerge from the presence of Joule heating and the temperature-dependence of the material properties and they are numerically solved by the Newton-Raphson algorithm. Notice that this material nonlinearity directly allows the inclusion of the Thomson effect. Furthermore, the formulation is dynamic and monolithic; the former feature is solved by backward finite differences and the latter is carried out by defining a coupled assembled matrix that increases the accuracy of the formulation. Finally, several examples are reported to show the capabilities of the numerical formulation that can be used as a “virtual laboratory.” In particular, h-refinements and residual convergence tests are conducted to validate the codification. Then, comparisons between analytical and numerical solutions for cooling thermoelectric cells are reported in order to highlight the advantages of the simulations against the simple one-dimensional analytical solutions. In conclusion, the use of the present numerical tool could be applied for a proper design and optimization of thermoelectric devices. For instance, this tool was used to optimize the shape of a pulsed thermoelectric in [25].