Thermoelectricity from Macro to Nanoscale: Wave Behaviour and Non-Local Effects

Abstract

In this chapter, the physical principles to be taken into account in thermoelectricity at the nanometre scale are discussed. We argue that the numerical methods must also be adapted to the emergence of new physical behaviours at that scale, namely, wave propagation of heat, diffusive-ballistic transition, nonlocal effects, among others. It is first shown that thermoelectric phenomena at the nanoscale can be described by introducing thermodynamic inertia and nonlocal effects. The transport equations are obtained from the thermodynamics of irreversible processes. After this, we introduce the Spectral Chebyshev Collocation method as a well-suited numerical method to deal with the new physical behaviours appearing at the nanoscale. We then show the use of these formalisms to analyse specific and interesting aspects of the optimization of pulsed thermoelectricity and coupled thermoelectric modules.

1. Introduction: irreversible thermodynamics of thermoelectricity

Heat, electron and hole transport and generation-recombination of electron-hole pairs are processes which determine the functioning of thermoelectric devices. In this section, the principles of irreversible thermodynamics of modelling of thermoelectric phenomena are exposed. The contributions to the entropy production in the stationary state due to the dissipative effects associated with electron and hole transport, generation-recombination of electron–hole pairs as well as heat transport are analysed.

The operation of thermoelectric devices occurs on the basis of the so-called cross effects. On the one hand, the Seebeck effect which arises when an external temperature difference is applied leading to a charge flux and, on the other hand, the Peltier effect which causes thermal fluxes in the presence of an applied voltage [1, 2]. This is also connected with the Joule and Thomson effects arising with the electrical current together with the non-equilibrium electrons and holes and electron–hole recombination phenomena. The latter becomes a thermal source causing an internal energy heterogeneity described through the internal energy balance equation [3].

The importance of electron–hole recombination must be remarked since the optimization of the device must be achieved in the regime of bipolar transport [4]. In this context, the problem of heat dissipation and its removal concerns the parameters that characterize the work of the device. When the dimension of the device goes to the nanometric length scale, nonlocal and memory effects must be taken into account [5]. This is achieved through the introduction of thermodynamic inertia in the constitutive equations of the dissipative fluxes and local dependence of transport coefficients, respectively. In the following, the constitutive equations of dissipative fluxes will be derived from the very principles of irreversible thermodynamics [6]. Then nonlocal and inertial effects will be introduced.

Let us consider the system constituted by two species (electron and holes) and the lattice. Firstly, define the electron and hole densities by $n\left(\overrightarrow{r},t\right)$ and $p\left(\overrightarrow{r},t\right)$ . The continuity equations for the electrons (electric charge $-q$ ) and holes (electric charge $+q$ ) are given by

$$q\frac{\partial n}{\partial t}-\nabla \cdot {\overrightarrow{J}}_n=-\mathit{qR},$$ E1

$$q\frac{\partial p}{\partial t}+\nabla \cdot {\overrightarrow{J}}_p=-\mathit{qR},$$ E2

where ${\overrightarrow{J}}_n$ and ${\overrightarrow{J}}_p$ are electric charge fluxes of electrons and holes, respectively, $q$ is the elementary electric charge and $R$ is the balance of generation and recombination of electron–hole pairs processes. Now we write the balance equation of total internal energy $u$ including the contribution of electrons $u_n$ , holes $u_p$ and the lattice $u_L$ :

$$\frac{\partial u}{\partial t}+\nabla \cdot {\overrightarrow{J}}^u=-\overrightarrow{E}\cdot \left({\overrightarrow{J}}_n-{\overrightarrow{J}}_p\right)+P,$$ E3

with $u=u_n+u_p+u_L$ and ${\overrightarrow{J}}^u={\overrightarrow{J}}_n^u+{\overrightarrow{J}}_p^u+{\overrightarrow{J}}_L^u$ the total internal energy flux. The total electric field (external plus self-consistent field) is represented by $\overrightarrow{E}$ . Finally, $P$ is a source term which can include light energy transference to the lattice and other processes. We use the expression for the Gibbs equation for each of the components of the system in order to find the balance equation of the total entropy density. We begin by writing the corresponding Gibbs equation for electrons, holes and the lattice. They are

$$\begin{array}{l}{T}_n{\mathit{ds}}_n={\mathit{du}}_n-{\mathrm{\Phi }}_n\mathit{dn},\\ T_p{\mathit{ds}}_p={\mathit{du}}_p+{\mathrm{\Phi }}_p\mathit{dp},\\ T_L{\mathit{ds}}_L={\mathit{du}}_L,\end{array}$$ E4

In these equations, ${\mathrm{\Phi }}_x$ is the electrochemical potential of species $x=n,p$. Explicit expressions for the electrochemical potentials are the following.

$$\begin{array}{l}{\mathrm{\Phi }}_n={\nu }_n-\mathit{q\phi },\\ {\mathrm{\Phi }}_p={\nu }_p+\mathit{q\phi }\end{array}$$ E5

being $\phi$ the total electric potential (external plus self-consistent field). At this point, it is convenient to mention that the self-consistent electric field is given by Poisson’s equation:

$$\nabla \cdot \left(\epsilon \nabla {\phi }_s\right)=q\left(N_A^{-}+N_D^{+}+n-p\right),$$ E6

where $\epsilon$ is the permittivity constant and $N_A^{-},N_D^{+}$ are the densities of ionized acceptors and donors, respectively. It is now assumed that the components of the system are in thermal equilibrium, that is, $T_n=T_p=T_L\equiv T$ . The sum of Eqs. (4) yields the balance equation for the total entropy density. One gets

$$\frac{\partial s_T}{\partial t}+\nabla \cdot {\overrightarrow{J}}_T^s={\overrightarrow{J}}_T^s\cdot \nabla \left(\frac{1}{T}\right)+\frac{1}{T}{J}_n\cdot {\overrightarrow{E}}_n^{\text{`},}+\frac{1}{T}{J}_p\cdot {\overrightarrow{E}}_p^{,}+R\left(\frac{{\mathrm{\Phi }}_n-{\mathrm{\Phi }}_p}{T}\right)+\frac{P}{T}.$$ E7

In obtaining Eq. (7) use has been made of Eqs. (1–3). The total entropy density $s_T$ in Eq. (7) is then given by.

$$s_T=s_n+s_p+s_L.$$ E8

The total entropy flux ${\overrightarrow{J}}_T^s$ has been defined as

$${\overrightarrow{J}}_T^s=\frac{1}{T}\left({\overrightarrow{J}}_q+q^{-1}{\mathrm{\Phi }}_n{\overrightarrow{J}}_n+q^{-1}{\mathrm{\Phi }}_p{\overrightarrow{J}}_p\right),$$ E9

and the heat flux ${\overrightarrow{J}}_q$ becomes:

$${\overrightarrow{J}}_q={\overrightarrow{J}}^u.$$ E10

The field ${\overrightarrow{E}}_x^{,}$ in Eq. (7) is a generalized electric field given by ${\overrightarrow{E}}_x^{,}=q^{-1}\nabla {\mathrm{\Phi }}_x\mp \overrightarrow{E}$ , with $x=n,p$ .

The right-hand side of Eq. (7), excepting the term $P/T$ , is the entropy generation term divided by the temperature, that is,

$$\mathit{T\sigma }={\overrightarrow{J}}_n\cdot {\overrightarrow{E}}_n^{\text{`},}+{\overrightarrow{J}}_p\cdot {\overrightarrow{E}}_p^{,}-{\overrightarrow{J}}_T^s\cdot \nabla T+R\left({\mathrm{\Phi }}_n-{\mathrm{\Phi }}_p\right),$$ E11

and it has the form $\sum {\overrightarrow{J}}_i\cdot {\overrightarrow{X}}_i$ , where ${\overrightarrow{J}}_i$ and ${\overrightarrow{X}}_i$ are generalized fluxes and thermodynamic forces, respectively. This allows us to identify the dissipative fluxes, and therefore, the origin of irreversibilities in the system. The generalized fluxes are $\left({\overrightarrow{J}}_n,{\overrightarrow{J}}_p,{\overrightarrow{J}}_T^s,R\right)$ and the corresponding forces $\left({\overrightarrow{E}}_n^{,},{\overrightarrow{E}}_p^{,},\nabla T,{\mathrm{\Phi }}_n-{\mathrm{\Phi }}_p\right)$ .

The second law of thermodynamics demands that $\sigma >0.$ This condition is satisfied if the thermodynamic forces and fluxes are linearly related as follows:

$$\left(\begin{array}{l}{{\displaystyle \overrightarrow{J}}}_n\\ {{\displaystyle \overrightarrow{J}}}_p\\ {{\displaystyle \overrightarrow{J}}}_T^s\\ R\end{array}\right)=\left(\begin{array}{llll}{L}_{\mathit{nn}}& 0& L_{\mathit{ns}}& 0\\ 0& L_{\mathit{pp}}& L_{\mathit{ps}}& 0\\ L_{\mathit{sn}}& L_{\mathit{sp}}& L_{\mathit{ss}}& 0\\ 0& 0& 0& L_{\mathit{RR}}\end{array}\right)\left(\begin{array}{l}{{\displaystyle \overrightarrow{E}}}_n^{,}\\ {{\displaystyle \overrightarrow{E}}}_p^{,}\\ \nabla T\\ {\mathrm{\Phi }}_n-{\mathrm{\Phi }}_p\end{array}\right),$$ E12

where the coefficients $L_{\mathit{xy}}$ are the so called Onsager coefficients. They are determined through phenomenological arguments and obey the reciprocity Onsager relations, namely, $L_{\mathit{xy}}=L_{\mathit{yx}}$ , with $x,y=n,p,s,R$ . Eq. (12) are the constitutive equations of the system which make complete the description offered by Eqs. (1–3, 6) together with the caloric equation $u_T=c_{V}T$ , where $c_V$ is the specific heat at constant volume. The constitutive equation, Eq. (12), contain well-known phenomenological laws: Ohm, Fourier, Fick, and Peltier and Seebeck effects. In the following section, we expose some additional considerations to be taken into account when the dimensions of the thermoelectric systems are in the nanometric length scale.

2. Non-local and memory effects

In this section, we address the problem of heat transport in a thermoelectric nanoscaled layer when an electric current circulates through it. At that length scale, nonlocal and memory effects must be included. The former are due to size effects on the transport coefficients when going to the nanometric scale and the second one become from the thermodynamic inertia of the system. Let us consider the thermoelectric system shown in Figure 1, and let us pay attention to one of the semiconductor branches, the $n$ type for instance.

Figure 1.

Thermoelectric couple system. From Ref. [5].

The analysis of heat transport in that element of the device of Figure 1 is based on the constitutive equations obtained in Section 1, Eqs. (12). The equations explicitly read [7, 8, 9]

$${\tau }_{\mathit{eff}}\frac{\partial q}{\partial t}+q=-\left[K\left(L\right)+\mathrm{\Pi }{S}_E\left(L\right)\sigma \left(L\right)\right]\frac{\partial T}{\partial x}-\mathrm{\Pi }\sigma \left(L\right)\frac{\partial V}{\partial x},$$ E13

$${\tau }_J\frac{\partial J}{\partial t}+J=-S_E\left(L\right)\sigma \left(L\right)\frac{\partial T}{\partial x}-\sigma \left(L\right)\frac{\partial V}{\partial x},$$ E14

where for the sake of simplicity, we have denoted the heat and electric charge fluxes as $q$ and $J$ , respectively. Inertia and no locality were introduced through the time derivatives of the heat and electric charge fluxes and the dependence of the transport and thermoelectric coefficients on the width of the layer, respectively. The equations have been written in one spatial coordinate since if the length $L$ goes to the nanometric scale, the dimensions of the branch $n$ in the $y$ and $z$ directions become much larger than $L$ . The times ${\tau }_{\mathit{eff}}$ and ${\tau }_J$ are the relaxation times of the heat and electric charge fluxes, respectively. $K$ , $\sigma$ and $S_E$ are the thermal conductivity, the electric conductivity and the Seebeck coefficient, respectively. In Eqs. (13) and (14), they are denoted to depend on the length $L$ . Expressions for the transport coefficients as functions of $L$ are obtained within the higher order dissipative fluxes framework of extended irreversible thermodynamics [10, 11]. The expressions are the following

$$K\left(L\right)=\frac{K_0{L}^2}{2{\pi }^2{l}_p^2}\left(\sqrt{1+4{\pi }^2{\left(\frac{l_p}{L}\right)}^2}-1\right),$$ E15

$$\sigma \left(L\right)=\frac{{\sigma }_0{L}^2}{2{\pi }^2{l}_e^2}\left(\sqrt{1+4{\pi }^2{\left(\frac{l_e}{L}\right)}^2}-1\right),$$ E16

$$S_E\left(L\right)=\frac{4{\pi }^2{S}_0{l}_e^2}{3L^2}{\left(\frac{2\pi \left(l_e/L\right)}{\text{arctan}\left(2\pi \left(l_e/L\right)\right)}-1\right)}^{-1},$$ E17

with $l_p$ and $l_e$ the mean free path of heat and electric charge carriers, respectively. $K_0$ , ${\sigma }_0$ and $S_0$ are the bulk thermal conductivity, electric conductivity and Seebeck coefficient, respectively. A hyperbolic-type transport equation for temperature can be obtained by introducing Eqs. (13)–(14) in Eq. (3) (with $P=0$ ) and using the caloric equation for the inner energy. The procedure may be followed in Appendix A of reference [5]. The resulting equation is

$${\alpha }_{\mathit{eff}}\frac{{\partial }^{2}T}{\partial t^2}+\frac{\partial T}{\partial t}=\nu \frac{{\partial }^{2}T}{\partial x^2}+\beta J^2+\mathit{\varsigma J}\frac{\partial }{\partial t}\left(\frac{\partial T}{x}\right),$$ E18

where the dimensionless coefficients are defined as follows:

$${\alpha }_{\mathit{eff}}=\frac{{\tau }_{\mathit{eff}}}{\tau },\nu =\frac{K\left(L\right)\tau }{\rho C_v{L}^2},\beta =\frac{\tau J_0^2}{\rho C_v\sigma \left(L\right)T_h},\varsigma =\frac{J_0{\tau }_{\mathit{eff}}{S}_E\left(L\right)}{\rho C_{v}L}.$$ E19

The characteristic time $\tau =\frac{{\pi }^2{L}^2\rho C_v}{4K\left(L\right)}$ is the diffusion time and $T_h$ is a reference temperature. Equation (18) is here solved numerically for a thermoelectric thin film (in the branch $n$ , see Figure 1) subjected to a Dirichlet boundary condition in the hot side of the thermoelectric device and a Robin type one in the cold side. We choose the space domain as $-1\le x\le 1$ , in such a way that boundary conditions can be written as:

$${\left.\frac{\partial T}{\partial x}\right|}_{x=-1}=\mathit{\gamma JT},T\left(1\right)=1,$$ E20

where the dimensionless coefficient $\gamma$ is defined as $\gamma =S_E\left(L\right)J_{0}L/\left(2K\left(L\right)\right)$ . As the initial condition, we state that the device is at room temperature, that is, $T\left(x,0\right)=1$ . Eq. (18) shows step-solutions which are very challenging numerically speaking. Thus, it was solved by using the numerical code based on the Spectral Chebyshev Collocation method described in Section 3. In Figure 2a, it is shown the time evolution of the cold temperature towards the stationary state. Since Silicon is a basic material for short and long scale devices, our departing results come from considering doped Silicon as working material, whose properties have been published before in [12]: $K_0=149$ Wm⁻¹ K⁻¹, ${\sigma }_0=35.5\times {10}^3$ Ω⁻¹ m⁻¹, $S_E=440\times {10}^{-6}$ V K⁻¹ and $a_E=88\times {10}^{-6}$ m² s⁻¹, where $a_E$ is the thermal diffusivity. The mean free path of the heat carriers and their mean velocity are [10] $l_p=40$ nm and $v=3K/\rho C_p{l}_p$ , respectively. The mean free path of electric charge carriers was assumed to be of the order of the lattice constant of Silicon $l_e=0.5$ nm. As hot side temperature, we use $T_h=373.1$  K [12]. We present the results without further discussion. The temperature starts to decrease once the electric current is applied. The transient to the stationary strongly depends on $L$ . As it can be appreciated in Figure 2a, as the Knudsen numbers $\delta$ (defined as $l_p/L$ ) increases the wave behaviour appears. On the contrary, for small $\delta$ , the steady state is reached quickly through a relaxation process without any oscillation. The minimum reached temperature is about $160K$ when $\delta =4$ . The response of the system to a short electric pulse superimposed to the stationary state obtained with the optimal electric current was also studied. The squared electric pulse with magnitude $3.5$ is applied at time $t=50$ . The wave behaviour of temperature produces subsequent supercooling transient process once the electric current pulse has been applied. The duration of it exceeds the duration of the supercooling obtained with the pulse. The minimum reached temperature is about 300 K when $\delta =4$ . In Figure 2a, it can also be seen an overheating after the application of the current pulse (which occurs at the dimensionless time 50). This effect is due to the thermal inertia, and it is more pronounced for small values of $L$ for which the thermal inertia plays a more important role. Figure 2b shows the evolution to the stationary state by omitting in the calculations the effect of the size $L$ on the Seebeck coefficient given by Eq. (17). The oscillating decay for small values of L remains but, interestingly, the minimum temperature reached is much higher than that reached when the effect of the size $L$ on the Seebeck coefficient is included by means of Eq. (17). Thus, the thermal efficiency of the cooling, measured by the temperature difference reached, is much greater due to the influence of the size on the Seebeck coefficient.

Figure 2.

(a) Cold side temperature $T_c$ as a function of time $t$ for different $\delta$ . From Ref. [5]. (b) Cold side temperature $T_c$ as a function of time $t$ for different L. In this case, there is no applied squared electric pulse and size effects on Seebeck coefficient have been omitted. From Ref. [13]. The initial temperature in both cases (a) and (b) is 373.1 K.

3. Computational methods

In this section, we present a couple of numerical methods in order to solve the hyperbolic transport equation (Eq. (18)) that models the heat transport in thermoelectric thin layers. Since Eq. (18) cannot be solved analytically, a numerical approach must be accomplished. Typically, low-order numerical methods such as finite differences (FD) are used. However, because of the high temporal gradients (due to the source term with crossed derivatives) for the small scales in the problem, that is, when the heat wave propagation or ballistic transport is present, the FD method fails in providing satisfactory results. As commented previously, this kind of equations is very challenging and we opt for a high-order numerical method in order to find its solution. High-order (or spectral) methods have previously been used to study heat transport based on the Maxwell-Cattaneo-Vernotte equation giving a hyperbolic transport equation in macroscopic systems [14, 15] and in some microscopic devices [16]. Thus, in this section, we discuss the numerical approximations of the one-dimensional hyperbolic heat transport equation for a low- and a high-order schemes.

3.1. Finite differences

The code for the finite differences scheme is a standard one. It considers a forward difference and a second-order central differences for the first and second derivatives in time, respectively. In turn, for the second-order spatial derivative, a second-order central differences was used as well. The grid points were uniformly spaced, the time integration was an explicit with constant time step.

3.2. Spectral Chebyshev collocation

The high-order numerical code is based on the Spectral Chebyshev Collocation (SCC) method. The method departs by establishing a partial sum of Chebyshev polynomials. The partial sum is then considered to represent the solution of a partial differential equation (PDE). The solution of the PDE equation is satisfied exactly at the Gauss-Lobatto collocation points

$$x_i=\text{cos}\left(\frac{\mathit{i\pi }}{N}\right),i=1,\dots ,N-1,$$ E21

where N denotes the number of points or the size of the grid. Thus, the solution is in the domain of the standard Chebyshev polynomials, that is {x| −1 < x < 1}. The partial sum of Chebyshev polynomials was inserted in the spatial derivatives of Eq. (18) obtaining expanded derivative matrices. The obtained coefficient equation system was solved by the matrix-diagonalization method in the physical space directly. A further explanation of the SCC method is found in [17, 18]. In order to compare directly with FD scheme, a coordinate transformation to interval to 0 < x < 1 was done. The time marching scheme was the same as the FD method.

3.3. Comparison

Figure 3a presents the steady-state temperature as a function of the spatial x-coordinate in the thermoelectric film for the micro-scale. We can easily observe that the temperature difference between the hot (T_h) and cold (T_c) sides is about 3.1 degrees. Because of the Joule effect, the temperature profile is parabolic. The cold side temperature T_c as a function of time is shown in Figure 4b. The departing point is the initial condition when T_c = 373.1 K, then the cooling Peltier effect acts till the temperature reaches the steady state T_c = 370 K. We can observe (in Figure 3) that both numerical schemes agree quantitatively in the steady state spatial distribution as well as modelling the transitory state. If we define the error between both solutions as

$$\epsilon =\text{max}\left|u_{\mathit{FD}}-u_{\mathit{SCC}}\right|.$$ E22

Figure 3.

Comparison of the numerical schemes in the micro-scale. (a) Steady-state temperature profile; (b) transient of the cold side. From Ref. [19].

Figure 4.

Comparison of the numerical schemes in the nanoscale. (a) Steady-state temperature profile; (b) transient of the cold side. From Ref. [19].

The maximum error for the steady profile (Figure 3a) with N = 30 is $\epsilon =8.6\times {10}^{-5}$ , and the error for the transitory is one order of magnitude higher, that is $\epsilon =2.5\times {10}^{-4}$ which is acceptable.

As the length of the system L diminishes, reaching the nanoscale, the steady state is a line with positive slope, see Figure 4a, whereas the wave heat transport is clearly visible during the transient, showing a damped harmonic oscillation, see Figure 4b. At this scale, the α_eff coefficient in the heat Eq. (18) becomes important, and thus, the wave term becomes dominant, see Table 1. The overall error for the steady-state solution is acceptable ( $\epsilon =9.8\times {10}^{-5}$ ). However, it is considerable larger ( $\epsilon =6\times {10}^{-3}$ ) for the transient. Finally, we can conclude that the SCC technique is more robust since the convergence of the solution is assured with smaller grid points than the FD method.

L(m)	${\alpha }_{\mathrm{eff}}$	$\beta$	$\gamma$
1 × 10^−4	2.16 × 10^−8	7.81 × 10^−3	7.38 × 10^−3
1 × 10^−5	2.16 × 10^−6	7.81 × 10^−5	7.38 × 10^−4
1 × 10^−6	2.16 × 10^−4	8.27 × 10^−7	7.82 × 10^−5
1 × 10^−7	2.16 × 10^−2	2.39 × 10^−8	2.26 × 10^−5
1 × 10^−8	2.16 × 10^0	2.00 × 10^−9	1.89 × 10^−5
1 × 10^−9	2.16 × 10^2	1.97 × 10^−10	1.85 × 10^−5

Table 1.

Dimensionless coefficients in Eqs. (18) and (20) as a function of the length of the system.

4. Optimal performance of thin thermoelectric layers

The control of heat in the nanoscale could have important consequences in applications as refrigeration, energy generation, energy transport and others. So, nanophonics has become a very active field of theoretical, computational and experimental research in the last 15 years. Many questions about the use of non-equilibrium thermodynamics principles at the micro and the nanoscale are being discussed, and several issues should be solved to make devices at those scales a matter of practical use. Particularly, devices in the micro and nanometric length scales work at high frequencies and they generate heat fluxes that can be in the order of thousands of watts per square centimetre. These irreversible processes elevate temperature reducing the device’s life time. This makes necessary the study of those operating conditions producing less dissipation, which often correspond to the minimum entropy production. Here we describe the effects of the width on the time evolution of temperature in thin thermoelectric layers and, particularly, on the thermal figure of merit and the entropy generation. The analysis is based on hyperbolic-type Eq. (18) describing the time evolution of dissipative flows including size effects on the thermal and electric conductivities. The hyperbolic Eq. (18) was solved with a numerical code explained in Section 3. The transition from the diffusive heat transport to the wave propagation regime is controlled by the system’s size when going from the micro to the nanometric scale of lengths.

When applying an electrical current to the thermoelectric system in order to obtain the maximum gradient between boundaries, the temperature distribution along the system’s length is parabolic as it can be seen in Figure 5a. Such current is named the optimal electrical current. In turn, the entropy generation is a parabolic decaying function showing its higher value at the cold end of the device. It is important to note that these profiles are invariant to the device length when applying the optimal current. Figure 5b shows the dependence of the temperature difference between the boundaries of the thermoelectric (Eq. (23)), the thermal figure of merit (Eq. (24)), the entropy generation (Eq. (25)), and, as functions of the length’s device (inset):

$$\mathit{\Delta T}=T_h-T_c,$$ E23

$${\mathit{ZT}}_h=\frac{2T_h\mathit{\Delta T}}{{\left(T_h-\mathit{\Delta T}\right)}^2},$$ E24

$$\underset{}{\overset{.}{S}}\left(x\right)=S_J\frac{J^2}{T}+S_F\frac{1}{T^2}{\left(\frac{\mathit{dT}}{\mathit{dx}}\right)}^2,$$ E25

where the coefficients are $S_J={J_0}^2/\sigma T_h$ and $S_F=K/l^2$ . In Figure 5b, we can note that devices with $L<{10}^{-6}$ m increase their thermal figure of merit, from ${\mathit{ZT}}_h=0.017$ (for $L<{10}^{-6}$ m) to ${\mathit{ZT}}_h\approx 0.4$ for $L<{10}^{-8}$ m.

Figure 5.

Temperature T and entropy generation S profiles. Lines and dots denote the exact and numerical solutions, respectively. $J_0=5.7785x{10}^7$ A/m2, $L=1.0x{10}^{-4}$ m. (b) Temperature difference between boundaries $\mathit{\Delta T}$ , thermal figure of merit ${\mathit{ZT}}_h$ and total entropy generation S as a function of L using the optimal electrical current. From Ref. [13].

We find that the steady-state thermal figure of merit (TFM) improves as the width goes to the nanoscale while the entropy generation increases. We identify a transition at a length of the order of the mean free path of heat carriers; in silicon it is about 40 nm. When going from the micro to the nanoscale, this transition is featured by an abrupt increment of both the total entropy production and the TFM. Above L = 40 nm, the heat transfer is dominated by diffusive processes; below this value, it is in the form of heat waves. The wave heat transport is clearly visible in the nonstationary process. An interesting study of the ballistic-diffusive transition in metals controlled by the wave number vector can be seen in [20].

5. Pulsed thermoelectric phenomena in thin films

As it was shown in Section 2, pulsed regimes produce a lower temperature than that obtained in the stationary state even with the optimal electric current for both uniform and non-uniform materials. This phenomenon is due to the fact that the Peltier effect occurs mainly at the cold junction while Joule heating is distributed in the bulk introducing a difference in the time taken by each one to influence the cold side of the device. The cold temperature is first changed by the Peltier effect and after diffusion Joule heat reaches the cold junction affecting it. Some examples of devices which need to be overcooled during a short time are mid-IR laser gas sensors [21], condensation hygrometers and microelectronic processors generating hotspots [22, 23, 24, 25, 26]. The effect of the pulse form has been widely studied in macroscale of lengths. It has been shown that by applying a quadratic pulse form, the supercooling effect can be improved over other forms [27]. Some other pulse forms present additional advantages [28, 29]. Here we explore the influence of the electric pulse shape in the supercooling effect when the dimensions of the thermoelectric device go to the submicrometre length scale. We study the effects of the shape of the electric pulse on the maximum diminishing of temperature by applying pulses in the form $t^a$ with a being a power going from 0 to 10.

In Figure 6, it can be seen the different shapes of the imposed electric pulse as a function of time. The duration in all cases is 0.163 and the maximum magnitude 3.5 over the stationary electric density (with normalized magnitude of one). The values are the optimal in order to obtain the maximum supercooling for the squared shape ( $t^0$ ). In the same figure, $t^0$ denotes the squared-shape pulse.

Figure 6.

Pulsed electric density current as a function of time. Different pulse shapes are used in optimizing the supercooling effect. From Ref. [30].

In Figure 7, it can be seen the time evolution of the temperature during the supercooling at the cold side $T_c$ of the thermoelectric for two distinct thickness of the film, namely, (a) $L=1\times {10}^{-4}$ m and (b) $L=1\times {10}^{-8}$ m. Each curve corresponds to one of the shaped pulses accordingly with the notation of Figure 6. The curves in Figure 7a reproduces the previous result found in [23] (microscopic case). The curves in Figure 7b are the result of the present analysis. It is remarkable the fact that the super cooling effect is about eight times larger at nanometric than at micrometric scale.

Figure 7.

Cold side temperature $T_c$ as a function of time t for due to pulse shapes: (a) Microscale and (b) nanoscale. Fractional pulse shapes perform better at the nanoscale. From Ref. [30].

Pulses with a fractionary number perform better for nanoscaled devices, whereas those with a bigger than unity do it for microscaled ones. We also find that the supercooling effect is improved by a factor of 6.6 over long length scale devices in the best performances and that the elapsed supercooling time for the nanoscaled devices equals the best of the microscaled ones.

6. Coupled thermoelectric modules

Thermoelectric systems are efficient devices for small size cooling objectives. Thus, this section is devoted to study the thermal performance of system composed of two thermoelectric devices. As a first step, a theoretical model of the heat transport in both thermoelectric devices will be developed. As a second step, the theoretical results will be compared with experimental data. Similar experimental devices have been previously reported in the literature [31, 32, 33, 34, 35, 36]. However, we use a local approach instead of global energy balances, which allow us to obtain spatial distributions of the main physical properties.

6.1. Experimental procedure

The experimental set-up is conformed by a two-stage Peltier cooler, that is, two thermoelectric modules connected electrically independent and thermally in series, see Figure 8. The modules (with side length L = 30 mm and width h = 3.6 mm) are denoted by M in the figure and are made of Bismuth Telluride alloys. The two-stage system is located on a metal plate (15 × 8 cm). The plate is kept at a constant temperature by the contact with hot water, which is continuously forced to circulate by a pump in a rectangular frame where the hot plate is located. The pump is a LMI MILTON ROY Microprocessor dosing and is denoted by B. The modules are operated by an electric current I from a BK PRECISION 1696 DC power supply V. The electric currents are in the range 0 < I < 1.517 A. Avoiding the transient, the temperature on top of module one $M_1$ was measured with a thermocouple (K type, Extech 470 True RMS Multimeter). Silitek thermal paste was used for the joints between the thermocouple and the cold wall, and the two modules. The latter is for avoiding the thermal decoupling of joints. When the two-stage system is turned on, the temperature in the cold wall is diminished because the modules generate a heat flow from the cold to hot wall. The error in the temperature data was obtained by adjusting a normal distribution to the data at same points.

Figure 8.

Scheme of two coupled thermoelectric devices. This study is devoted to the thermal analysis of the thermoelectric materials denoted by n type. The devices are actioned by injecting an electric current $I_i$ . The thickness (width) is denoted by $L_i$ . The i denotes the first (1) and second (2) device. From Ref. [37].

6.2. Mathematical model

The heat transport problem can be reduced by considering only the heat flux in a branch (n type) of the thin thermoelectric modules, see Figure 9. An electric current I_i, i = 1, 2, is injected through each of the thermoelectric modules. The system is thus subjected to a Dirichlet boundary condition in the hot side at $T_h$ , and a Robin type one in the cold side $T_c$ , see Eq. (20).

Figure 9.

Scheme of two coupled thermoelectric devices. This study is devoted to the thermal analysis of the thermoelectric materials denoted by n type. The devices are actioned by injecting an electric current $I_i$ . The thickness (width) is denoted by $L_i$ . The i denotes the first (1) and second (2) device. From Ref. [37].

At this point we assume that each of the thermoelectric modules behaves independent, as the single thermoelectric case (Eq. (18)), and they only share a common boundary. The heat transport for the coupled system is modelled by a system of two one-dimensional differential equations, each for every thermoelectric module

$$0={\nu }_1\frac{d^2{T}_1}{{\mathit{dx}}^2}+{\beta }_1{J_1}^2,0={\nu }_2\frac{d^2{T}_2}{{\mathit{dx}}^2}+{\beta }_2{J_2}^2.$$ E26

The subscript denotes the first (1) and second (2) modules. The boundary conditions for the cold and hot sides are, respectively,

$${\left[\frac{{\mathit{dT}}_1}{\mathit{dx}}\right]}_{x=-1}={\gamma }_1{J}_1{T}_1\left(-1\right),T_2\left(1\right)=1.$$ E27

While the coupling boundary conditions are

$${\gamma }_1{J}_1{T}_1\left(0\right)-{\left[\frac{{\mathit{dT}}_1}{\mathit{dx}}\right]}_{x=0}=\mathit{\text{δLδK}}{\gamma }_2{J}_2{T}_2\left(0\right)-\mathit{\text{δLδK}}{\left[\frac{{\mathit{dT}}_2}{\mathit{dx}}\right]}_{x=0},T_1\left(0\right)=T_2\left(0\right),$$ E28

being $\mathit{\delta L}=L_1/L_2$ and $\mathit{\delta K}=K_2/K_1$ the ratios between the lengths and the thermal conductivities of modules 1 and 2, respectively. The coupling boundaries at x = 0 are obtained from a heat balance and by equalling the temperatures of both devices. The solution to the system of Eqs. (26) reads

$$T_1\left(x\right)=a_1{x}^2+b_{1}x+c_1,T_2\left(x\right)=a_2{x}^2+b_{2}x+c_{2,}$$ E29

where the coefficients $a_i$ , $b_i$ and $c_i$ are given by.

$$\begin{array}{c}{a}_i=-\frac{{\beta }_i{J_i}^2}{2{\nu }_i},c_1=\frac{\mathit{\text{δLδK}}\left(a_2-1\right)+\frac{a_1\left(2+{\gamma }_1{J}_1\right)}{1+{\gamma }_1{J}_1}}{\frac{{\left({\gamma }_1{J}_1\right)}^2}{1+{\gamma }_1{J}_1}-\mathit{\text{δLδK}}\left(1+{\gamma }_2{J}_2\right)}=c_2,\\ b_2=1-a_2-c_2,b_1=c_1\left({\gamma }_1{J}_1-\mathit{\text{δLδK}}{\gamma }_2{J}_2\right)+\mathit{\delta L\delta }{\mathit{Kb}}_2.\end{array}$$ E30

The difference of temperature between the hot and cold side is

$$\Delta T_{\mathit{\text{coupled}}}=T_2\left(1\right)-T_1\left(-1\right).$$ E31

6.3. Comparison

Before presenting and discussing the theoretical results, the experimental measurements for a single and a coupled system are shown in Figure 10a, b. Figure 10a shows the temperature difference between the hot and cold sides $\mathit{\Delta T}$ as a function of the circulating electric current for a single thermoelectric device. A parabolic behaviour with a maximum $\mathit{\Delta T}\approx 25$  K can be appreciated. Figure 10b presents elliptic isocurves of $\mathit{\Delta T}$ as function of the two electric currents, in the case that the system is composed of two thermoelectric modules. For this array, the maximum of the temperature difference is $\mathit{\Delta T}\approx 34$  K, that is, nine degrees of extra cooling are obtained when using two thermoelectric modules instead of a single one. For both cases, the single and two thermoelectric module system, the electric current is normalized with $I_0=0.93$  A, which is the optimal current for a single module.

Figure 10.

Temperature gradient $\mathit{\Delta T}$ as function of the electric current. (a) and (c) Single device. (b) and (d) Coupled devices. First row, experimental measurements from reference [37]. Second row, theoretical predictions from Ref. [37].

Figure 10c shows the temperature difference $\mathit{\Delta T}$ between the hot and cold sides as a function of the electric current for a single thermoelectric device. This result can be compared with the experiment, Figure 10a, and the parabolic behaviour with a maximum is found. Figure 10d makes evident the dependence of $\mathit{\Delta T}$ when the system is composed on two thermoelectric devices, presenting elliptic isocurves. Theoretical results (Figure 10c, d) obtained from the simple model discussed before agree qualitatively with the experimental ones (Figure 10a, b) showing a maximum $\mathit{\Delta T}$ for certain values of the electric currents. Since the theoretical model used to get Figure 10c, d is a one-dimensional representation of a three-dimensional problem, a quantitatively comparison is far from reflecting the experimental measurements. However, if we calculate the percentage of the maximum extra cooling obtained when using two thermoelectric modules instead of a single one, that is $\left(\Delta T_{\mathit{\text{coupled}}}-\Delta T_{\text{sin}\mathit{gle}}\right)/\Delta T_{\text{sin}\mathit{gle}}$ , we find that it is the same for both the experimental and the theoretical results. The maximum extra cooling that could be reached is 36%.

Some conclusions derived from the above are the following. First, the theoretical solutions of the temperature difference between the hot and the cold sides show a good qualitative agreement with the experimental measurements. Second, an optimal cooling with respect to the electric currents circulating through the modules has been found. Finally, an improvement of 36% in the performance, measured as mentioned above, of the coupled thermoelectric modules with respect to a single thermoelectric module is theoretically predicted. The analysis made in this section can be useful for the design of thermoelectric coolers.

7. Discussion and conclusions

We have analysed different aspects of the performance of thermoelectric films when an electric current flows through them, maintaining one of their sides at a constant temperature. First, we considered the effects of the size of the film on the Seebeck coefficient and the thermal and electrical conductivities as well as thermal inertia. This forms the basis for all the analysis of the problem of reducing the dimensions to the nanoscale. The introduction of such effects was done through the use of Eqs. (15)–(17) in the constitutive equations of the irreversible thermodynamics of electric charge and heat transport, Eqs. (12). When the latter are combined with the conservation equations of charge and energy, Eqs. (1)–(3), we arrive at the heat transport equation that was used systematically throughout the development of the analysis, Eq. (18). The different effects studied, namely, the wave behaviour of the propagation of heat, the response of the material to a pulse of electric current, the coupling of thermoelectric couples, the transition from the diffusive regime of heat transport to wave propagation, and so on were explained in terms of Eq. (18). The conclusions were diverse: (1) thermodynamic inertia improves the thermal performance of a cooler, (2) the transition from diffusive transport to wave propagation of heat is controlled by the size of the thermoelectric material, (3) in the operating pulsed mode in a cooler, the shape of the electrical pulse is crucial for its thermal performance, (4) it is always possible to find values ​​of the thermal parameters of a device that optimize its performance, (5) the operation of thermoelectric systems is closely related to the production of entropy in the system. Behind all these behaviours is the combined effect of the reduction of the thermal and electrical conductivities of the material and the increase of the Seebeck coefficient with the reduction of the dimensions. To close this section we show in Figure 11 the behaviour of the thermoelectric figure of merit, $\mathit{ZT}=\sigma S_E^{2}T/K$ , when the thickness of a doped silicon semiconductor film decreases towards the nanoscale of lengths. As can be seen in the figure, the increase in the Seebeck coefficient dominates the reduction of thermal and electrical conductivities. The net effect is the increase of the thermoelectric figure of merit. Notably, the figure at a length of $l={10}^{-9}m$ is greater than that at $l={10}^{-5}m$ by a factor of $287$ .

Figure 11.

Effect of Eqs. (15)–(17), which describe the variation of the transport and thermoelectric coefficients with the size of the material $l$ , on the thermoelectric figure of merit defined as usual: $\mathit{ZT}=\sigma S_E^{2}T/K$ .

8. Perspectives

Undoubtedly, a work perspective in thermoelectric power generation and refrigeration systems is to bring to the nanoscale the dimensions of the elements that make up structured systems. There is a lot of work done in systems composed or structured from elements in the millimetric escale of lengths or in larger scales [38], but a very promising perspective is the construction of those energy converting systems from nano-components such as for instance, the case of 1D phononic crystals. Several problems on thermoelectric phenomena that require attention when going from the macro to the nanoscale are: effects of size on the thermoelectric properties of the components of the nanostructured systems, the wave character of the propagation of heat and the phenomena of thermal interference and resonance, the effect of wave propagation of heat in pulsed systems, the effects of size on irreversible processes, the simultaneous consideration of the effects of size and the temperature dependence of the properties of materials, the effect of thermal inertia on the wave propagation of heat. It is also necessary to study the relationship between entropy generation and thermoelectric conversion efficiency and between structure and efficiency. In relation to this, a working hypothesis is that it is possible to find a structure that minimizes the production of entropy and to exploit interference and resonance phenomena to improve the thermoelectric energy conversion efficiency.

Acknowledgments

F.V. and A.F. acknowledge financial support from CONACYT-México under contracts 183358 and 258623, respectively. F.V. also thanks the project “Fronteras de la Ciencia” 367. A. F. also thanks the Cátedras program from CONACYT.

Conflict of interest

We declare that there is no conflict of interest.