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


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

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

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.

• 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

the interpolation functions.

274 Bringing Thermoelectricity into Reality

3. Governing equations

Historical note 1

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 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 Ω, boundary Γ with outward normal n and its surrounding Ω<sup>∞</sup>. A general expression of a balance equation for the quantity ϒ at time t and in any material point x is given by:

$$\frac{d}{dt}\Big|\_{\Omega}\rho(\mathbf{Y};\mathbf{x},t)\quad d\varOmega\quad=-\int\_{\Omega}\nabla\cdot\mathbf{J}(\mathbf{Y};\mathbf{x},t)\quad d\varOmega\quad+\int\_{\Omega}\Xi(\mathbf{Y};\mathbf{x},t)\quad d\varOmega\quad,\tag{1}$$

where r, J and Ξ denote the ϒ-density, ϒ-current density and ϒ-production, respectively. This equation balances the quantity ϒ since it contains information on:


In thermoelectricity, both electric and thermal fields are present and, consequently, there are two ϒ quantities as shown in Table 2. As observed, the quantities are the free electric charge r f q and the specific entropy s; the current densities are the electric j and entropy j<sup>s</sup> fluxes. Regarding the production terms, there is not production for the electric field due to the fact that the


Table 2. Quantities to be balanced in thermoelectricity.

electric charge can neither be created nor destroyed. Consequently, the electric balance becomes a conservation equation. On the contrary, the entropy production Σ<sup>s</sup> is due to the irreversibilities (also called dissipations) generated by both the Fourier heat transport and the Joule heating.

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 \dot{\jmath}, \qquad \rho\_m \dot{\mathbf{s}} = -\nabla \cdot \dot{\jmath}\_s - \frac{\dot{\jmath} \cdot \nabla V}{T}, \tag{2}$$

j ¼ �γð Þ T ∇V � αð Þ T γð Þ T ∇T, q ¼ �κð Þ T ∇T þ αð Þ T Tj: (6)

Computational Thermoelectricity Applied to Cooling Devices

http://dx.doi.org/10.5772/intechopen.75473

277

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,

In semiconductors, the material properties γ, κ and α commonly depend on temperature. For instance, for the metal-metalloid alloy Bi2Te3 the dependency on temperature of these properties is shown in Figure 2. As observed, both γ and κ exhibit a strong temperature dependency that results in a material nonlinearity. Explicitly, the temperature dependency of these material

, <sup>κ</sup>ð Þ¼ <sup>T</sup> <sup>κ</sup><sup>0</sup> <sup>þ</sup> <sup>κ</sup>1<sup>T</sup> <sup>þ</sup> <sup>κ</sup>2T<sup>2</sup>

Notice that the temperature-dependency of αð Þ T is the responsible of the Thomson's effect and, consequently, the inclusion of material nonlinearities is a requirement for a proper modeling of

The boundary conditions are additional constrains at the boundary Γ that are present at all the boundary value problems. These conditions guarantee the existence of a unique solution and,

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

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

therefore, the problem is well posed. There exist two types of boundary conditions:

, <sup>α</sup>ð Þ¼ <sup>T</sup> <sup>α</sup><sup>0</sup> <sup>þ</sup> <sup>α</sup>1<sup>T</sup> <sup>þ</sup> <sup>α</sup>2T<sup>2</sup>

: (7)

based on the degree Celsius, was also other of his achievements.

properties can be fitted to second order polynomials to give:

<sup>γ</sup>ð Þ¼ <sup>T</sup> <sup>γ</sup><sup>0</sup> <sup>þ</sup> <sup>γ</sup>1<sup>T</sup> <sup>þ</sup> <sup>γ</sup>2T<sup>2</sup>

thermoelectricity.

3.3. Boundary conditions

freedom at the boundary.

the fluxes at the boundary.

Historical note 2

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 j<sup>s</sup> ¼ q=T, where q is the heat flux to give:

$$\rho\_m \dot{\mathbf{s}} = -\nabla \cdot \left(\frac{\mathbf{q}}{T}\right) \underbrace{-\frac{\mathbf{q} \cdot \nabla T}{T^2} - \frac{\dot{\mathbf{j}} \cdot \nabla V}{T}}\_{\Sigma\_s}.\tag{3}$$

Now, the classical heat balance equations can be obtained by manipulating (3) and by taking into account that <sup>s</sup>\_ <sup>¼</sup> cT\_ <sup>=</sup>T, where <sup>c</sup> is the specific heat capacity, to give:

$$
\rho\_m \circledast \dot{T} = -\nabla \cdot \boldsymbol{q} - \dot{\boldsymbol{j}} \cdot \nabla V. \tag{4}
$$

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 ∇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 Σ<sup>s</sup> 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:

$$\begin{Bmatrix} -\frac{\nabla V}{T} \\\\ -\frac{\nabla T}{T^2} \end{Bmatrix} = \begin{bmatrix} \frac{\alpha^2}{\kappa} + \frac{1}{\gamma T} & -\frac{\alpha}{\kappa T} \\\\ -\frac{\alpha}{\kappa T} & \frac{1}{\kappa T^2} \end{bmatrix} \begin{Bmatrix} j \\\\ q \end{Bmatrix} \tag{5}$$

where γ, κ and α denote electric and thermal conductivities and Seebeck coefficient, respectively. Finally, the set of two coupled equations of thermoelectricity becomes:

$$\dot{q} = -\gamma(T)\nabla V - a(T)\gamma(T)\nabla T,\qquad \mathbf{q} = -\kappa(T)\nabla T + a(T)T\mathbf{j}.\tag{6}$$

#### Historical note 2

electric charge can neither be created nor destroyed. Consequently, the electric balance becomes a conservation equation. On the contrary, the entropy production Σ<sup>s</sup> is due to the irreversibilities (also called dissipations) generated by both the Fourier heat transport and the Joule heating.

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 equa-

<sup>0</sup> <sup>¼</sup> <sup>∇</sup> � <sup>j</sup>, <sup>r</sup>ms\_ ¼ �<sup>∇</sup> � <sup>j</sup><sup>s</sup> � <sup>j</sup> � <sup>∇</sup><sup>V</sup>

where the second term on the right-hand side of the right equation represents the Joule heating

The entropy balance of (2) can be manipulated by considering that j<sup>s</sup> ¼ q=T, where q is the heat

Now, the classical heat balance equations can be obtained by manipulating (3) and by taking

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 ∇V

Within the nonequilibrium thermodynamic formalism [31], the transport equations are obtained from the entropy production Σ<sup>s</sup> 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:

> α2 κ þ 1 <sup>γ</sup><sup>T</sup> � <sup>α</sup> κT

� α κT

where γ, κ and α denote electric and thermal conductivities and Seebeck coefficient, respec-

1 κT<sup>2</sup>

8 >>>>><

>>>>>:

j

9 >>>>>=

>>>>>;

q

� <sup>q</sup> � <sup>∇</sup><sup>T</sup>

<sup>T</sup><sup>2</sup> � <sup>j</sup> � <sup>∇</sup><sup>V</sup> T

<sup>r</sup><sup>m</sup> <sup>c</sup> <sup>T</sup>\_ ¼ �<sup>∇</sup> � <sup>q</sup> � <sup>j</sup> � <sup>∇</sup>V: (4)


T � �

<sup>r</sup>ms\_ ¼ �<sup>∇</sup> � <sup>q</sup>

into account that <sup>s</sup>\_ <sup>¼</sup> cT\_ <sup>=</sup>T, where <sup>c</sup> is the specific heat capacity, to give:

� <sup>∇</sup><sup>V</sup> T

9 >>>=

>>>; ¼

tively. Finally, the set of two coupled equations of thermoelectricity becomes:

8 >>><

>>>:

� <sup>∇</sup><sup>T</sup> T2

<sup>T</sup> , (2)

: (3)

, (5)

tions in local forms read:

276 Bringing Thermoelectricity into Reality

term, as can be observed from (6).

3.2. Transport equations

flux to give:

and V, T denote voltage and temperature, respectively.

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 γ, κ and α commonly depend on temperature. For instance, for the metal-metalloid alloy Bi2Te3 the dependency on temperature of these properties is shown in Figure 2. As observed, both γ and κ 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:

$$\boldsymbol{\gamma}(T) = \boldsymbol{\gamma}\_0 + \boldsymbol{\gamma}\_1 T + \boldsymbol{\gamma}\_2 T^2, \quad \boldsymbol{\kappa}(T) = \boldsymbol{\kappa}\_0 + \boldsymbol{\kappa}\_1 T + \boldsymbol{\kappa}\_2 T^2, \quad \boldsymbol{a}(T) = a\_0 + a\_1 T + a\_2 T^2. \tag{7}$$

Notice that the temperature-dependency of αð Þ T 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 Γ 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:


Figure 2. Temperature-dependency of material properties for the Bi2Te3 material. Picture taken from [21].

In thermoelectricity, these conditions read:

$$\begin{aligned} \text{Dirichlet type}: & \quad V = \overline{V}, \quad & T = \overline{T}, \\ \text{Neumann type}: & \quad j \cdot \mathbf{n} = \overline{j}, \quad q \cdot \mathbf{n} = \overline{q}. \end{aligned} \tag{8}$$

i) The continuum domain Ω is divided into subdomains called finite elements Ω<sup>e</sup> interconnected at nodal points a, see Figure 3. Mathematically, this discretization is expressed by Ω ≈ ∪nel

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-

ii) The values at nodes are called degrees of freedom g and are the problem unknowns. In

<sup>g</sup><sup>a</sup> <sup>¼</sup> <sup>V</sup><sup>~</sup> <sup>a</sup>; <sup>T</sup>~a; ~\_

iii) A set of functions denominated "shape functions" 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

<sup>x</sup> <sup>≈</sup> <sup>N</sup>ax~a, V <sup>≈</sup> <sup>N</sup> <sup>a</sup> <sup>V</sup><sup>~</sup> a, T <sup>≈</sup> <sup>N</sup> aT~a, <sup>T</sup>\_ <sup>≈</sup> <sup>N</sup> <sup>a</sup> ~\_

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

a) The balance equations of (2) (left) and (4) are multiplied by test functions of the degrees of

Ta

n o: (9)

Computational Thermoelectricity Applied to Cooling Devices

http://dx.doi.org/10.5772/intechopen.75473

Ta, (10)

dimensional eight noded elements, as that shown in Figure 3, are considered.

thermoelectricity, they are voltage V, temperature T and their time derivatives:

the numerical interpolations become:

freedom δV and δT.

where the Einstein summation convention is used.

probably the most systematic approach; it consists of three steps:

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

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

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

<sup>e</sup>¼<sup>1</sup>Ωei . 279

where V, T, j and 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 Ω is divided into subdomains called finite elements Ω<sup>e</sup> interconnected at nodal points a, see Figure 3. Mathematically, this discretization is expressed by Ω ≈ ∪nel <sup>e</sup>¼<sup>1</sup>Ωei . 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, threedimensional eight noded elements, as that shown in Figure 3, are considered.

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

$$\mathbf{g}\_a = \left\{ \tilde{V}\_a, \tilde{T}\_a, \dot{\tilde{T}}\_a \right\}. \tag{9}$$

iii) A set of functions denominated "shape functions" 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:

$$\mathbf{x} \approx \mathcal{N}\_a \tilde{\mathbf{x}}\_{a\nu} \qquad V \approx \mathcal{N}\_a \, \tilde{V}\_{a\nu} \qquad T \approx \mathcal{N}\_a \tilde{T}\_{a\nu} \qquad \dot{T} \approx \mathcal{N}\_a \dot{\tilde{T}}\_{a\nu} \tag{10}$$

where the Einstein summation convention is used.

In thermoelectricity, these conditions read:

278 Bringing Thermoelectricity into Reality

4. Finite element formulation

using the following steps:

Historical note 3

mechanics.

respectively.

Dirichlet type : V ¼ V, T ¼ T, Neumann type : j � n ¼ j, q � n ¼ q,

Figure 2. Temperature-dependency of material properties for the Bi2Te3 material. Picture taken from [21].

where V, T, j and q denote prescribed voltage, temperature, electric flux and thermal flux,

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

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

(8)

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 δV and δT.


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

$$\begin{aligned} \begin{bmatrix} \int\_{\Omega} \nabla \delta V \cdot & \mathbf{j} d\Omega \ -\oint\_{\Gamma} \delta V \overline{\mathbf{j}} \ \mathrm{d}\Gamma &=& 0, \\ \int\_{\Omega} [\nabla \delta T \cdot \mathbf{q} + \delta T (\mathbf{j} \cdot \nabla V + \rho\_m c \, \dot{T})] & \, d\Omega \ -\oint\_{\Gamma} \delta T \overline{\mathbf{q}} \ \mathrm{d}\Gamma &=& 0. \end{aligned} \tag{11}$$

Finally, the solution is updated by using g<sup>k</sup>þ<sup>1</sup>

to the University of California at Berkeley (USA).

extracted, voltage and temperature distributions are analyzed.

where Tc and Th refer to the temperature at cold and hot ends, respectively.

5. Numerical examples

to half of the thermocouple.

<sup>b</sup> <sup>¼</sup> <sup>g</sup><sup>k</sup>

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

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

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 Qc <sup>¼</sup> <sup>82</sup>:01 (W) with voltage drop <sup>V</sup> = 15.33 (V) and under Th <sup>¼</sup> Tc <sup>¼</sup> <sup>50</sup><sup>∘</sup>

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

Figure 4. Left: Sketch of the thermoelectric cooling module. Right: Finite element mesh and boundary conditions applied

<sup>b</sup> <sup>þ</sup> <sup>d</sup>g<sup>k</sup>

<sup>b</sup> until the convergence is reached. In

http://dx.doi.org/10.5772/intechopen.75473

Computational Thermoelectricity Applied to Cooling Devices

C,

281

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

$$\begin{aligned} \delta V & \approx \mathcal{N}\_a \delta \tilde{V}\_a & \delta T \approx \mathcal{N}\_a \delta \tilde{T}\_a \\ \nabla \delta V & \approx \nabla \mathcal{N}\_a \delta \tilde{V}\_a = \mathcal{B}\_a \delta \tilde{V}\_a & \nabla \delta T \approx \nabla \mathcal{N}\_a \delta \tilde{T}\_a = \mathcal{B}\_a \delta \tilde{T}\_a. \end{aligned} \tag{12}$$

Introducing these expressions in the weak form of (11), the FE residuals for each finite element of domain Ω<sup>e</sup> of boundary Γ<sup>e</sup> become:

$$\begin{array}{rcl} \mathcal{R}\_{b}^{V} &=& \int\_{\Omega\_{\epsilon}} \mathcal{B}\_{b}^{\mathsf{T}} \boldsymbol{\jmath} \, \mathrm{d}\Omega\_{\epsilon} - \oint\_{\Gamma\_{\epsilon}} \mathcal{N}\_{b} \, \boldsymbol{\jmath} \, \mathrm{d}\Gamma\_{\epsilon\prime} \\ \mathcal{R}\_{b}^{\mathsf{T}} &=& \int\_{\Omega\_{\epsilon}} \left[ \mathcal{B}\_{b}^{\mathsf{T}} + \mathcal{N}\_{b} \left( \boldsymbol{\jmath} \cdot \boldsymbol{\nabla} \boldsymbol{V} + \rho\_{m} \boldsymbol{c} \, \boldsymbol{\bar{T}} \right) \right] \, \, \mathrm{d}\Omega\_{\epsilon} - \oint\_{\Gamma\_{\epsilon}} \mathcal{N}\_{b} \, \boldsymbol{\eta} \, \, \mathrm{d}\Gamma\_{\epsilon\prime} \end{array} \tag{13}$$

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 Δ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}\_{\mathfrak{a}}\right|^{k} = -\frac{\left|\partial\mathcal{R}\_{\mathfrak{a}}\right|^{k}}{\left|\mathfrak{d}\mathfrak{g}\_{\mathfrak{b}}\right|}\,^{k}\mathrm{dg}\_{\mathfrak{b}}^{k} \Rightarrow -\frac{\left|\partial\mathcal{R}\_{\mathfrak{a}}\right|^{k}}{\left|\mathfrak{d}\mathfrak{g}\_{\mathfrak{b}}\right|}\,^{k} = \left|\mathcal{L}\_{1}\mathcal{K}\_{\mathfrak{a}\mathfrak{b}} + \mathfrak{c}\_{2}\,\mathcal{L}\_{\mathfrak{a}\mathfrak{b}\mathfrak{b}}\right|\tag{14}$$

where a, b refer the local numbering of two generic FE nodes, the parameters c<sup>1</sup> and c<sup>2</sup> depend on the time integration scheme, see [19], and the consistent tangent matrices K and the nonzero capacity matrix C are derived at each iteration as:

$$\begin{split} \mathcal{K}\_{ab}^{VV} &= -\frac{\partial \mathcal{R}\_a^V}{\partial \dot{V}\_b} = -\int\_{\Omega\_c} \mathcal{B}\_a^T \frac{\partial \dot{f}}{\partial \dot{V}\_b} \quad \text{d}\Omega\_{\epsilon\prime} \\ \mathcal{K}\_{ab}^{VV} &= -\frac{\partial \mathcal{R}\_a^V}{\partial \dot{T}\_b} = -\int\_{\Omega\_a} \mathcal{B}\_a^T \frac{\partial \dot{f}}{\partial \dot{T}\_b} \quad \text{d}\Omega\_{\epsilon\prime} \\ \mathcal{K}\_{ab}^{TV} &= -\frac{\partial \mathcal{R}\_a^T}{\partial \dot{V}\_b} = -\int\_{\Omega\_a} \left[ \mathcal{B}\_a^T \frac{\partial \mathcal{q}}{\partial \dot{V}\_b} + \mathcal{N}\_a \frac{\partial (\dot{f} \cdot \nabla V)}{\partial \dot{V}\_b} \right] \quad \text{d}\Omega\_{\epsilon\prime} \\ \mathcal{K}\_{ab}^{TT} &= -\frac{\partial \mathcal{R}\_a^T}{\partial \dot{T}\_b} = -\int\_{\Omega\_a} \left[ \mathcal{B}\_a^T \frac{\partial \mathcal{q}}{\partial \dot{T}\_b} + \mathcal{N}\_a \frac{\partial (\dot{f} \cdot \nabla V)}{\partial \dot{V}\_b} \right] \quad \text{d}\Omega\_{\epsilon\prime} \\ \mathcal{C}\_{ab}^{TT} &= -\frac{\partial \mathcal{R}\_a^T}{\partial \dot{T}\_b} = -\int\_{\Omega\_c} \mathcal{N}\_a \, \rho\_m \, \text{c} \; \mathcal{N}\_b \quad \text{d}\Omega\_{\epsilon\prime} \end{split} \tag{15}$$

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

Finally, the solution is updated by using g<sup>k</sup>þ<sup>1</sup> <sup>b</sup> <sup>¼</sup> <sup>g</sup><sup>k</sup> <sup>b</sup> <sup>þ</sup> <sup>d</sup>g<sup>k</sup> <sup>b</sup> 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).
