# An Analytical Development of the Hyperbolic Behaviour of Micro Thermoelectric Coolers

Authors: Giulia De Aloysio; Giampaolo D’Alessandro; Filippo de Monte
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

The transient behaviour of a micro thermoelectric cooler is described in the present paper through the theory of the thermal wave, involving the relaxation time. The internal heat generation due to the sudden application of the electric current is taken into account by means of the Heaviside function. The governing equations of both the semiconductors are solved by a modified Separation of Variables Method that allows us to have a better description of the device which can be obtained at early times. As regards the performances, the cooling load shows discontinuities due to the contributions of the travelling wave fronts of both the semiconductors. The results show that the coefficient of performance (COP) of the device obtained with the hyperbolic model is lower than that provided by the parabolic model at early times.


---

## Body

## 1. Introduction

In the literature many models related to the transient behaviour of thermoelectric devices have been presented [1–9]. These systems have large application as coolers [9–12], thanks to the possibility of maintaining the junction temperature as low as required. Most of the models make use of the parabolic thermal diffusion equation: three-dimensional numerical coupled models have been proposed recently, in which the coupling of the thermal and the electric field is considered, [4, 10]. In particular, the coupling of both the fields is analyzed under unsteady and steady states for coolers in [4, 10].However, the parabolic modelling is inadequate for devices with small dimensions, in particular in the analysis of their transient behaviour, during the start-up phase and the shut-down phase and when changing the operational parameters. In all of these cases it is more adequate to use the hyperbolic heat conduction model, which considers the relaxation time of the material [6–8]. This is a crucial parameter in the heat conduction theory because it takes into account the finite speed of the heat propagation. The relaxation time is absent in the classical theory of heat diffusion based on the Fourier law, whose intrinsic underlying assumption is that heat propagates in the material with an infinite speed.Unlike the previously proposed model, in the present paper we describe the sudden application of the electric current by means of the Heaviside function. Following Haji-Sheikh and Beck’s approach [13], the thermal field in both the semiconductors is obtained analytically by means of a modified Separation of Variables (SoV, for short) Method for time-independent boundary conditions of the first kind. According to this method, the solution of the hyperbolic problem is obtained by modifying the solution of the corresponding parabolic problem through an appropriate time-dependent function. This function is due to the presence of the relaxation time and may be derived by imposing the modified solution which will be the solution of the hyperbolic problem. Then MATLAB ambient is used to implement the solution. Once the temperature distributions are known, the performances of the device, cooling load, and COP are evaluated. The results of the analysis show that the COP obtained by means of the hyperbolic model is lower than that provided by the parabolic model at early times.


## 2. Analytical Formulation

The elementary unit of a thermoelectric cooler (TEC) module consists of two semiconductors, P and N, connected electrically in series, as shown in Figure1. Because of the Peltier effect, by applying an electric current, heat is absorbed at the cold junction, whose temperature is denoted by Tc, and heat is released at the hot junction, whose temperature is Th. The heat absorption at the cold junction is counteracted by Joule and Fourier’s effects. It is worth noting that the Thomson effect is neglected in this analysis. The following assumptions are made to perform the analytical transient modelling of the device:(a)
The thermal conductivity and the electric resistivity are considered temperature-independent and computed at an average temperatureTavg=(Th+Tc)/2, by using (23) (see ahead). This assumption is valid for low electric applied currents, as shown in [10].(b)
Adiabatic side surfaces of the TEC are taken into account.(c)
Semiconductors are employed in the analysis, which have the same geometrical characteristics, the cross-sectional areasA1 and A2, and the lengths of the thermoelements, L1 and L2. We assume that L1=L2=L.(d)
At the instantt=0 both the semiconductors have an ambient uniform temperature T0 and a sudden constant uniform current is applied to the system.Figure 1
1D schema of a thermoelectric cooler.


### 2.1. Governing Equations for the Semiconductors

The energy equations for both the semiconductors are derived as suggested in [6]. However, in the present paper the energy equation is modified, as instantaneous internal heat generation has been considered by means of the Heaviside function in the hyperbolic model. Therefore, the energy balance equation valid for each semiconductor may be written as (1)ρc∂T∂t+∂q∗∂x+eJ∂Φ∂xHt=0.Ht represents the Heaviside function. In particular it is(2)Ht=0ift<01ift>0.The thermal flux at time t+τ- is defined as in (3a) as follows:(3a)q∗x,t+τ-=-k∂T∂x∓TSeJHt,where τ- is the relaxation time and the minus sign stands for semiconductor N, whereas the plus sign stands for semiconductor P.By applying Taylor expansion to (3a) and by neglecting the higher order terms, (3a) may be written as(3b)q∗x,t+τ-∂q∗∂t=-k∂T∂x∓TSeJHt.By combining (1) and (3b) and by performing some other calculations, as suggested in [6], the governing equation that allows us to determine the temperature distribution under the hyperbolic heat conduction model for both the pellets is given by (4a), (4b), and (4c). Consider (4a)kA∂2T∂x2+RI2LHt-ρAc∂T∂t-τ-ρAc∂2T∂t2+τ-∂∂tRI2LHt∓SI∂T∂xHt=0.Performing the derivative of the last term in (4a), a mixed derivative term should be generated. This term is neglected because its contribution is significantly lower than the one due to the other terms in (4a), if the length scale of the pellet is of the order of millimeters, as suggested in [8]. In order to solve (4a), initial and boundary conditions are provided in(4b)Tx,0=T0,∂T∂tx,0=0,(4c)T0,t=Tc,TL,t=Th.Before solving the problem, it is convenient to write it in dimensionless form as in [6], by introducing some dimensionless parameters, listed in the Nomenclature. Therefore, the dimensionless governing equation for both the semiconductors and the initial and the boundary conditions may be expressed by (5a), (5b) and (5c). Consider (5a)Fi∂2ϑi∂ξ2+Di-∂ϑi∂η-τi∂2ϑi∂η2+τiDiδη=0,(5b)ϑiξ,0=0,∂ϑi∂ηξ,0=0,(5c)ϑi0,η=1,ϑi1,η=Th-T0Tc-T0=h.The subscripti indicates the type of the semiconductor; in fact i=1 is used for semiconductor N, whereas i=2 is used for semiconductor P. The dimensionless groups Di and Fi are defined as(6a)D1=LR1I2A1k1Tc-T0,D2=L2R2I2CRA2k1L2Tc-T0,(6b)F1=1,F2=αR.The parabolic heat diffusion equations are easily obtained by considering that for the classical heat diffusion equation the relaxation time disappears in (5a); that is, τi=0 and only one initial condition on the temperature is needed to determine the temperature distribution inside the pellets.


## 2.1. Governing Equations for the Semiconductors

The energy equations for both the semiconductors are derived as suggested in [6]. However, in the present paper the energy equation is modified, as instantaneous internal heat generation has been considered by means of the Heaviside function in the hyperbolic model. Therefore, the energy balance equation valid for each semiconductor may be written as (1)ρc∂T∂t+∂q∗∂x+eJ∂Φ∂xHt=0.Ht represents the Heaviside function. In particular it is(2)Ht=0ift<01ift>0.The thermal flux at time t+τ- is defined as in (3a) as follows:(3a)q∗x,t+τ-=-k∂T∂x∓TSeJHt,where τ- is the relaxation time and the minus sign stands for semiconductor N, whereas the plus sign stands for semiconductor P.By applying Taylor expansion to (3a) and by neglecting the higher order terms, (3a) may be written as(3b)q∗x,t+τ-∂q∗∂t=-k∂T∂x∓TSeJHt.By combining (1) and (3b) and by performing some other calculations, as suggested in [6], the governing equation that allows us to determine the temperature distribution under the hyperbolic heat conduction model for both the pellets is given by (4a), (4b), and (4c). Consider (4a)kA∂2T∂x2+RI2LHt-ρAc∂T∂t-τ-ρAc∂2T∂t2+τ-∂∂tRI2LHt∓SI∂T∂xHt=0.Performing the derivative of the last term in (4a), a mixed derivative term should be generated. This term is neglected because its contribution is significantly lower than the one due to the other terms in (4a), if the length scale of the pellet is of the order of millimeters, as suggested in [8]. In order to solve (4a), initial and boundary conditions are provided in(4b)Tx,0=T0,∂T∂tx,0=0,(4c)T0,t=Tc,TL,t=Th.Before solving the problem, it is convenient to write it in dimensionless form as in [6], by introducing some dimensionless parameters, listed in the Nomenclature. Therefore, the dimensionless governing equation for both the semiconductors and the initial and the boundary conditions may be expressed by (5a), (5b) and (5c). Consider (5a)Fi∂2ϑi∂ξ2+Di-∂ϑi∂η-τi∂2ϑi∂η2+τiDiδη=0,(5b)ϑiξ,0=0,∂ϑi∂ηξ,0=0,(5c)ϑi0,η=1,ϑi1,η=Th-T0Tc-T0=h.The subscripti indicates the type of the semiconductor; in fact i=1 is used for semiconductor N, whereas i=2 is used for semiconductor P. The dimensionless groups Di and Fi are defined as(6a)D1=LR1I2A1k1Tc-T0,D2=L2R2I2CRA2k1L2Tc-T0,(6b)F1=1,F2=αR.The parabolic heat diffusion equations are easily obtained by considering that for the classical heat diffusion equation the relaxation time disappears in (5a); that is, τi=0 and only one initial condition on the temperature is needed to determine the temperature distribution inside the pellets.


## 3. Thermal Field: Solution Method

The analytical solution of (5a) may be obtained by using the SoV technique as in [13].The hyperbolic temperature field for both the pellets is provided by (7a) and (8a). It is expressed as the sum of a steady state part ϑss,iξ and of a complementary part:(7a)ϑ1ξ,η=ϑss,1ξ+exp⁡-η2τ1∑n=1∞An,1sin⁡λnξΨn,1η+exp⁡-η2τ1∑n=1∞CnD1sin⁡λnξζn,1η,(7b)ϑss,1ξ=-D12ξ2+h+D12-1ξ+1,(7c)ζn,1η=sinh⁡βn,12-λn4ηβn,12-λn4,(8a)ϑ2ξ,η=ϑss,2ξ+exp⁡-η2τ2∑n=1∞An,2sin⁡λnξΨn,2η+exp⁡-η2τ2∑n=1∞CnD2sin⁡λnξζn,2η,(8b)ϑss,2ξ=-D22αRξ2+h+D22αR-1ξ+1,(8c)ζn,2η=sinh⁡βn,22-λn4αR2ηβn,12-λn4αR2.The eigenvaluesλn=nπ are the same for both the semiconductors. Moreover,(9a)βn,i=λn2Fi-12τi,(9b)Ψn,iη=cosh⁡βn,i2-λn4Fi2η+ζn,i2τi,(9c)An,i=2λn-1nh+Di-1n-1λn2Fi-1,(9d)Cn=21--1nλn.The parabolic solution for both the semiconductors may be obtained as suggested in [14] and is expressed by(10a)ϑ1ξ,η=ϑss,1ξ+∑n=1∞An,1sin⁡λnξexp⁡-λn2η,(10b)ϑ2ξ,η=ϑss,2ξ+∑n=1∞An,2sin⁡λnξexp⁡-αRλn2η.The subscripti always refers to the kind of semiconductor: i=1 for semiconductor N and i=2 for the P one.


## 4. Heat Fluxes: Solution Method

In order to complete the thermal analysis of the transient behaviour of the device, after obtaining the thermal fields, it is necessary to determine the conductive heat fluxes of each semiconductor. As regards the hyperbolic model, the thermal fluxes are given according to the non-Fourier heat equation, as in [15]. In fact the conductive heat flux for each pellet may be expressed as in (11a)qix,t+τi-∂qix,t∂t=-kiAi∂Tix,t∂xwith the initial condition(11b)qix,0=0.By means of the dimensionless variables, listed in the Nomenclature, (11a) and (11b) may be rewritten as (12a) and (12b). Consider (12a)Qiξ,η+τi∂Qiξ,η∂η=-vi∂ϑiξ,η∂ξ,(12b)Qiξ,0=0,where v1=1 for semiconductor N and v2=kRAR for semiconductor P. By integrating (12a) and by applying the initial condition, the hyperbolic thermal fluxes for each semiconductor are provided by(13)Qiξ,η=-viτiexp⁡-ητi∫0η∂ϑiξ,ω∂ξξexp⁡ωτidω.As regards the parabolic heat fluxes, the relaxation time disappears in (11a) and each dimensionless conductive contribution may be expressed as in(14)Qiξ,η=-vi∂ϑiξ,η∂ξ.After evaluating the temperature gradient, according to (7a), it is substituted in (13). By integrating it and after some other calculations, the dimensionless hyperbolic conductive heat flux for each semiconductor is provided by (15) and (16). Consider(15)Q1ξ,η=h+D112-ξ-1exp⁡-ητ1-1-exp⁡-η/2τ1τ1∑n=1∞Bn,1cos⁡λnξζn,1η+∑n=1∞En,1,(16)Q2ξ,η=v2h+D2αR12-ξ-1exp⁡-ητ2-1-v2exp⁡-η/2τ2τ2∑n=1∞Bn,2cos⁡λnξζn,2η+v2∑n=1∞En,2,(17)Bn,i=An,iλn+CnDi2λnFi,(18)En,i=CnDicos⁡λnξλnFiexp⁡-η2τicosh⁡βn,i2-λn4Fi2η-CnDicos⁡λnξλnFiexp⁡-ητi.


## 5. Performances of the Microcooler

The cooling and the heating load effects are expressed by the thermal fluxesQc and Qh, respectively. They may be written in dimensionless terms by considering an energy balance for the hot and the cold junctions, as suggested by(19)Qc=G+Q1,c0,η+Q2,c0,η,(20)Qh=M+Q1,h1,η+Q2,h1,η.The dimensionless termsQi,j represent the conductive heat fluxes for semiconductors N and P, provided, in the hyperbolic case, by (15) and (16) and evaluated at the cold and the hot ends, ξ=0 and ξ=1, respectively. G and M represent the dimensionless contributions due to the Peltier effect and they are defined in the Nomenclature. In particular, G is the thermal power absorbed at the cold junction and M is the thermal power released at the hot junction.In the parabolic model (21a) and (21b) has to be considered as follows: (21a)Qc=G-∂ϑ10,η∂ξ-kRAR∂ϑ20,η∂ξ,(21b)Qh=M-∂ϑ11,η∂ξ-kRAR∂ϑ2LR,η∂ξ.It can be seen that the cooling load and the thermal power released at the hot junction, provided by (19)–(21b), are valid for a single thermoelectric pair. In the case of n-pairs, it is sufficient to multiply the above mentioned equations for the number of pairs. The performances of a thermoelectric cooler are evaluated by means of the coefficient of performance, COP, defined in(22)COP=QcQh-Qc.


## 6. Results and Discussion

A miniature thermoelectric cooler as in [4, 10] has been considered for the simulations. The geometrical characteristics and some properties of both the semiconductors are listed in Table 1. As regards the thermoelectric materials, Bi2(Te0.94Se0.06)3 and (Bi0.25Sb0.75)Te3 are selected for the N- and P-type semiconductors, respectively, as in [4, 10]. The temperature Tc of the cold end of the semiconductors is set at 27°C, whereas the hot ends temperature is set at Th=67°C as in [10], where the temperatures are given in Kelvin. The working current is assumed to be 1 A. The initial temperature is set at 20°C as in [10].Table 1
Geometrical sizes and properties of the semiconductors.
Parameter
N-type
P-type
Reference
L (m)
1 × 10−3
1 × 10−3
[4, 10]
A (m2)
0.25 × 10−6
0.25 × 10−6
[4, 10]
ρ (kg/m3)
10 922.08
10 922.08
[3]
c (J/kg K)
160
190
[3]In order to consider the thermal conductivity, the Seebeck coefficient, and the electrical resistivity of the above mentioned materials, a quadratic relationship, described by (23), is employed. Consider(23)pj=pjTref1+ajT-Tref+bjT-Tref2,where pj indicates the generic property, that is, k, S, and ρe, which is evaluated, in the present paper at an average temperature Tavg=(Th+Tc)/2. Tref is provided in [10] and the coefficients aj and bj are listed in [4, 10]. The subscript j=1 is used for the thermal conductivity, j=2 is used for the Seebeck coefficient, and j=3 is used for the electrical resistivity. The relaxation time value considered for the simulations of the hyperbolic thermal behaviour is 20 ms.Figures2 and 3 show the temperature profiles in semiconductor N for different time values. The wave fronts travel in the material at a speed of α/τ-.Figure 2
Temperature profiles for semiconductor N for early time values (a). Enlargement of the wave front (b). Relaxation time = 20 ms.
(a)
(b)Figure 3
Temperature profiles for semiconductor N (high time values). Relaxation time = 20 ms.During their motion, these wave fronts give constructive interference and then they are reflected by the opposite end of the semiconductor.By the observation of the plotted temperature profiles, it is evident that for early time values the temperature profiles have oscillations at the wave fronts.This phenomenon is due to the numerical instability caused by the impossibility of representing the infinite terms in the series in (7a) and (8a). To remove this numerical instability, the Separation of Variables Method, which has been used here for solving the hyperbolic equation and coming from [13], should be replaced by a different approach. For example, we might use a finite series solution at early times coming from the combination of Laplace transform and method of images (e.g., sources and sinks) as given in [16]. But this is a subject for future research.It is worth noting that the above mentioned oscillations are more evident in semiconductor P, Figures4 and 5, due to the different diffusivity of this pellet, more precisely, because of the presence of the thermal diffusivity ratio in (8a).Figure 4
Temperature profiles for semiconductor P (early time values). Relaxation time = 20 ms.
(a)
(b)Figure 5
Temperature profiles for semiconductor P (high time values). Relaxation time = 20 ms.For time values higher than 230 ms, the hyperbolic temperature field becomes similar to the parabolic one and the oscillations are not so marked as in the early times as depicted in Figures3 and 5.Figure6 shows the influence of the relaxation time on the temperature profiles for semiconductor N. It is clear that if the relaxation time decreases, the temperature profile becomes similar to the parabolic one, as the speed of the thermal wave in the semiconductor increases. Moreover, at the thermal wave regions, a larger relaxation time value affects the term τiDiδη significantly in (5a). The conductive heat flux in semiconductor N, provided by (15) and evaluated at the cold end ξ=0, is represented in Figure 7(a). Some discontinuity points are evident in the plot of the hyperbolic contribution whose enlargement is given by Figure 7(b).Figure 6
Influence of the relaxation time on the temperature profile for semiconductor N whent=60 ms.Figure 7
Conductive heat flux at the cold end of semiconductor N (a). Enlargement of the discontinuities (b).
(a)
(b)This phenomenon takes place because each reflection of the wave front produces discontinuities in the heat flux. For early time values, the parabolic model does not allow the heat flux at the cold end to be properly represented. In fact, it goes to infinity when the time approaches zero. For time values higher than 500 ms the hyperbolic and the parabolic heat fluxes tend to overlap.In Figure8 the conductive heat flux at the hot end of semiconductor N has been depicted.Figure 8
Conductive heat flux at the hot end of semiconductor N (a). Enlargement of the discontinuities (b).
(a)
(b)Figure9 shows the cooling load profile, whose analytical expression is provided by (19). In comparison with Figure 7, the cooling load shows further discontinuities due to the contributions of the wave fronts due to the second semiconductor (the P one).Figure 9
Cooling load (a). Enlargement of the discontinuities (b).
(a)
(b)This occurs because the two wave fronts travel at different speeds, due to the different thermal diffusivities. It has to be observed that the value of the cooling load att=0 is not zero, because of the contribution of the Peltier effect.As regards the COP, it is plotted for time values higher than 200 ms. In fact, for earlier time values, the employed model, based on imposed time-independent boundary temperatures, does not allow us to describe accurately the real thermal behaviour of the device: the thermal power released at the hot junction and the cooling load cross as depicted in Figure10, providing a zero denominator in (22).Figure 10
Thermal power released at the hot junction and cooling load. Relaxation time 20 ms.As the time increases, the hyperbolic COP is similar to the parabolic one, reaching the value of about 0.72 as shown in Figure11.Figure 11
COP as a function of time.Figure12 shows the COP profile as a function of the time for different values of the parameter F2 (thermal diffusivity ratio). By increasing this ratio, it can be observed that the COP increases at large time values. In Figure 13 the COP behaviour as a function of the time is plotted for different values of the parameter D1. This parameter expresses the dimensionless internal heat generation due to the Joule effect. By increasing D1, the volumetric effects, due to the heat dissipation, become more and more relevant, if compared with the Peltier heat, which is a superficial effect. This phenomenon causes a sensible reduction of the cooling load, which affects the COP. The parameter D2 has the same influence on the COP, as the one depicted in Figure 13.Figure 12
COP as a function of time for different values ofF2.Figure 13
COP as a function of time for different values ofD1.Moreover, as regards the COP, the values deriving from the present analysis under steady state conditions are in agreement with those presented in the numerical parabolic coupled model available in the literature [10], as shown in Table 2. However, by applying an electric current of 1.5 A, it can be observed that the COP value under steady state condition is noticeably higher than the value presented in the numerical dynamic parabolic coupled model [4]. This is due to the more realistic time-dependent boundary conditions employed in the simulations.Table 2
Comparison of the COP values under steady state conditions.
I (A)
Present
Reference [10]
Reference [4]
0.6
1
0.98
—
0.9
0.80
0.80
—
1.2
0.56
0.60
—
1.5
0.45
0.43
0.22It is worth noting that in the coupled numerical parabolic model [4] the steady state is reached after about 600 s due to the higher thermal inertia of the system. Therefore, the future commitment is focussed on the hyperbolic modelling developed by means of the more realistic time-dependent boundary conditions of the first kind. This is done by considering the presence of the heat sink at the hot end of the system.


## 7. Conclusion

An analytical development of the hyperbolic behaviour of a micro thermoelectric cooler has been presented in which the sudden application of the electric current is considered by means of the Heaviside function. The thermal fields, obtained by means of the Separation of Variables technique with imposed boundary conditions of the first kind, show oscillations by the wave fronts for very early times, which are more marked for semiconductor P. The reflection of the wave fronts of both the semiconductors yields discontinuities in the cooling load profiles.Moreover, the not realistic imposed time-independent boundary conditions of the first kind allow us to represent the COP profile only for time values higher than 0.2 s, as the intersection of the thermal power released at the hot junction and the cooling load occurs.The future commitment is focussed on the hyperbolic modelling developed by means of the more realistic time-dependent boundary conditions of the first kind. This is done by considering the presence of the heat sink at the hot end of the system.



---

*Source: 697639-2015-11-08.xml*