A three-dimensional theoretical model for predicting transient thermal behavior of thermoelectric coolers

International Journal of Heat and Mass Transfer 53 (2010) 2001–2011

Contents lists available at ScienceDirect

International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt

A three-dimensional theoretical model for predicting transient thermal behavior of thermoelectric coolers Chin-Hsiang Cheng a,*, Shu-Yu Huang a, Tsung-Chieh Cheng b a b

Department of Aeronautics and Astronautics, National Cheng Kung University, Tainan 70101, Taiwan, ROC Department of Mechanical Engineering, National Kaohsiung University of Applied Sciences, Kaohsiung 807, Taiwan, ROC

a r t i c l e

i n f o

Article history: Received 14 August 2009 Received in revised form 3 December 2009 Accepted 3 December 2009 Available online 6 February 2010 Keywords: Thermoelectric cooler Theoretical model Transient behavior Experiment

a b s t r a c t A simulation model is developed and used to predict transient thermal behavior of the thermoelectric coolers. The present model amends the previous models, in which the P–N pair is simply treated as a single bulk material so that the temperature difference between the semiconductor elements was not possible to evaluate. Based on the present simulation model, the thermoelectric cooler is divided into four major regions, namely, cold end (region 1), hot end (region 2), and the P-type and N-type thermoelectric elements (regions 3 and 4). Solutions for the three-dimensional temperature fields in the P-type and the N-type semiconductor elements and transient temperature variations in the cold and the hot ends have been carried out. The magnitude of the coefficient of performance (COP) of the thermoelectric cooler are calculated in wide ranges of physical and geometrical parameters. To verify the numerical predictions, experiments have been conducted to measure the temperature variations of both the cold and the hot ends. Close agreement between the numerical and the experimental data of the temperature variations has been observed. Ó 2009 Elsevier Ltd. All rights reserved.

1. Introduction

eral, a thermoelectric cooler consists of a number of semiconductor element pairs which are connected electrically in series and thermally in parallel, and each pair includes a P-type and an N-type elements. The P-type and the N-type elements are used to ensure that the carriers flow in the same direction. Although in theory a single piece of semiconductor material could be adopted, the series connection is used to avoid the high electrical current requirement with the single element. An applied current causes charged carriers in the material, whether they are electrons or holes, to diffuse from the cold side to the hot side. These mobile charged carriers migrating to the hot side leave behind their oppositely charged and immobile nuclei at the cold side, thus giving rise to a temperature difference. Since a temperature difference creates a driving potential for the phonons to drift back to the cold side resulting in heat conduction at equilibrium. Improving the performance of a thermoelectric materials involves controlling the motion of phonons and electrons/holes. The most important performance index of the thermoelectric materials is the figure of merit (ZT) which is defined by

Thermoelectric effects in semiconductors cause electrical currents to flow due to temperature gradients and also cause temperature gradients when an electrical current is applied. A feature of the thermoelectric devices is the absence of moving parts when they are used to convert heat into electrical energy, and vice versa. Therefore, they can work over a long period of time in high performance without extra maintenance. This could be a significant advantage from system simplification and reliability points of view. Besides, these devices produce no waste matter in the conversion process. Therefore, it is regarded as an attractive method for energy conversion. Thermoelectric devices are basically categorized into two operation modes based on the direction of energy conversion: thermoelectric cooler (TEC) converting electricity to thermal energy and thermoelectric generator (TEG) converting heat into electricity. The thermoelectric cooler is a device in which an electric current is applied to semiconductor devices to produce an appreciable temperature difference at the two ends of the semiconductor. The colder end may then be utilized for cooling purposes. In gen-

ZT ¼ a2 T r=k

* Corresponding author. Address: Department of Aeronautics and Astronautics, National Cheng Kung University, No.1, Ta Shieh Road, Tainan 70101, Taiwan, ROC. Tel.: +886 6 2757575x63627. E-mail address: [email protected] (C.-H. Cheng).

where a is the Seebeck coefficient; T is the temperature; r is the electrical conductivity; k is thermal conductivity of the material [1,2]. A higher ZT value in TECs leads to higher COP, and in TEGs it leads to higher thermal efficiency. For the purpose of reaching higher figure of merit, one should try to elevate the magnitude of

0017-9310/$ - see front matter Ó 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.ijheatmasstransfer.2009.12.056

ð1Þ

2002

C.-H. Cheng et al. / International Journal of Heat and Mass Transfer 53 (2010) 2001–2011

Nomenclature A Ah As C h ha I k L La Lb M Q q Qc Qh

cross-sectional area of P-type or N-type element (m2) surface area of the hot end (cm2) total heat transfer surface area of the heat sink (m2) specific heat (J kg1 K1) heat transfer coefficient on surface of the heat sink (W m2 K1) heat transfer coefficient on side surfaces of P-type or N-type element (W m2 K1) electrical current (A) thermal conductivity (W m1 K1) height of the thermoelectric elements (m) length of the cross-sectional areas of P-type or N-type element in y-direction (m) length of the cross-sectional areas of P-type or N-type element in z-direction (m) mass (kg) cooling load (W) cooling load per unit area (W cm2), Q/Ah heat transfer rate between cold end and P-type or N-type element (W) heat transfer rate between hot end and P-type or N-type element (W)

the electrical conductivity and reduce the thermal conductivity [3]. Besides, a higher electrical conductivity of the material decreases the electrical resistance for electrical currents and hence lowers the Joule heating, and a lower thermal conductivity of the material increases the thermal resistance for heat conduction from the hot end to the cold end. To the authors’ knowledge, Stilbans and Regel [4] are the earliest researchers who carried out the dynamic behavior of the thermoelectric cooler. They found that the cooling rate is decreased as the thickness of the thermoelectric elements is increased. Parrot [5] extended the study to include the effects of variable heat conductivity, which is treated as a function of temperature. Gray [6] derived a mathematical model for predicting the dynamic behavior of TEC; however, the author did not take the hot end’s dynamic behavior into consideration. Bywaters and Blum [7] proposed the equivalent stage method (ESM) which was derived by simplifying the governing equations of Gray’s model, and used it to estimate the dynamic behavior of a multi-layer TEC. Huang and Duang [8,9] presented a theoretical model, utilized the numerical solutions of the model to control the temperature difference of the TEC device, and introduced the influences of thermal resistance into their model for calculating the efficiency of TEC. Recently, a dimensionless mathematical model is constructed by Naji et al. [10] to investigate the influences of the dimensionless parameters on the performance of TEC. Authors also found that the cold end temperature was elevated as the electrical resistivity is decreased and was reduced as the Thomson coefficient is increased. Xuan et al. [11] constructed a general model to study the effects of electrical resistances in the internal and external interface layers on the performance of the thermoelectric devices. Huang et al. [12] addressed the effects of the Thomson effect, the Joule heating, the Fourier’s heat conduction, and the radiation and convection heat transfer on the temperature distribution. Authors also found that the cooling efficiency of the thermoelectric cooler can be improved not only by increasing the figure of merit (ZT) but also by taking advantage of the Thomson effect. Khire et al. [13] combined the thermoelectric devices with the solar energy units to design an active building envelope (ABE) system and mentioned that the thermal resistance of the heat sink is essential to determine the

T Tc Th T1 V x, y, z

temperature (K) temperature of the cold end (K) temperature of the hot end (K) ambient temperature (K) volume of the thermoelectric elements (m3) rectangular coordinates

Greek symbols Seebeck coefficient of thermoelectric elements (V K1) electrical resistivity of thermoelectric elements (X m) Thomson coefficient of thermoelectric elements (V K1) density of thermoelectric elements (kg m3)

a e l q

Subscripts b cold plate of cold end cp ceramic plate of cold end hp ceramic plate of hot end n N-type thermoelectric element p P-type thermoelectric element s heat sink

number of TEC units required. Chakraborty et al. [14] presented a thermodynamic modeling of a solid state thermoelectric cooling device based on the temperature-entropy analysis. Most recently, Lee and Kim [15] carried out the numerical analysis to evaluate the coefficient of performance of the thermoelectric micro-cooler with a nearly three-dimensional model which is used to investigate the effects of the temperature difference and the current. According to the literature survey [16], it is recognized that mostly the previous theoretical models are limited to the onedimensional problems. However, in these existing models, the P–N element pair is simply treated as a single bulk material so that the difference in thermal behavior between the two semiconductor elements was not possible to evaluate. Only a few studies were conducted for the three-dimensional ones, for example, the model of Lee and Kim [15]. In practices, the P-type and the N-type single crystal silicon elements are made by doping boron and phosphorous, respectively, into the silicon wafer. According to the existing information, the properties of the P-type and the N-type silicon are dependent on the material doped, dopant density and doping time. For example, Beadle et al. [17] addressed that the electrical resistivity of N-type silicon is dependent on the dopant density. When the N-type silicon is doped at a dopant density of 5  1017 cm3, its electrical resistivity can reach 2.3  104 X m. The properties of the two semiconductor elements can be rather different under different doping conditions. As a result, it may lead to significant difference in temperature response in the P-type and N-type elements. Under these circumstances, a three-dimensional model is developed in this study that treats the P-type and the N-type semiconductor elements as two separate parts and considers the difference in temperature between the two elements. The model is also used to predict transient thermal behavior of the thermoelectric coolers.

2. Theoretical model The schematic of a thermoelectric cooler containing one semiconductor element pair is shown in Fig. 1. A thermoelectric cooler consists of a number of semiconductor element pairs. Neverthe-

C.-H. Cheng et al. / International Journal of Heat and Mass Transfer 53 (2010) 2001–2011

Ambient Temperature T∞

hAs(Th-T ∞ )

Heat sink Region 2 (Hot-end, Th) Qhp

IαpTh

P

Regions 3 and 4

Ceramic plate Qhn

-Iα nTh

N Ceramic plate Qcp

Qcn

IαnTc

IαpTc

Region 1 (Cold-end, Tc)

Cold plate

P

lected and compared with the numerical predictions to verify the theoretical model. The theoretical model is developed based on the following assumptions: (1) The hot and cold ends of the thermoelectric cooler are assumed to be lumped-capacity systems inside which the temperature is changed with time but is uniformly distributed. (2) The material of the P-type and N-type elements is homogeneous and isotropic. Furthermore, all the thermal and electrical properties of the elements are constant and independent of temperature. (3) The heat transfer coefficients on the outer surfaces are constant, and the ambient temperature (T1) is fixed at 298 K. The governing equations for the present study are expressed in the following statement:

Q (a) Energy balance for all regions

2.1. Cold end (region 1) As shown in Fig. 1(a), the cold end consists of a ceramic plate and a cold plate. The cooling load is the heat transfer rate that is absorbed into the cold plate from the external object that is cooled by this TEC. By treating the ceramic plate and the cold plate as a system, the energy equation is expressed as

ðM cp C cp þ M b C b Þ

N

L

L xp

xn

zp yp

Lap

Lbp

2003

zn yn

Lbn

Lan

(b) Dimensions of P-type and N-type thermoelectric elements Fig. 1. Schematic of a thermoelectric cooler.

less, the thermal characteristics among the pairs are periodic, and hence, only one pair is considered here. In this study, the solution domain is divided into four major regions, namely, cold end (region 1), hot end (region 2), and the P-type and N-type thermoelectric elements (regions 3 and 4), as shown in Fig. 1(a). Numerical solutions for the three-dimensional temperature fields in the P-type and the N-type semiconductor elements and transient temperature variations in the cold and the hot ends can be carried out. In addition, based on the obtained temperature data, the magnitude of the coefficient of performance (COP) of the thermoelectric cooler are calculated in wide ranges of physical and geometrical parameters, including electrical current (I), cooling load per unit area (q), heat transfer coefficients on the surfaces (h and ha) and geometrical parameters such as total area of heat sink (As) and height of semiconductor elements (L). The P-type and the N-type elements are rectangular columns of equal volume. A heat sink is attached on the hot end of the thermoelectric cooler for heat dissipation. Note that in this study the temperature distributions within the P-type and N-type elements are predicted by the three-dimensional model. This is the major concern of the present paper. However, the hot and cold ends of the thermoelectric cooler are minor parts of the modules and their specifications are directly dependent on the practical applications. Therefore, they are assumed to be lumpedcapacity bodies to simplify the mathematical model. For validation of the solution model, experiments are also conducted for a certain particular cases. Measurement data regarding the temperature variations in the cold and the hot ends are col-

dT c ¼ Q þ Q cp þ Q cn  Iðap  an ÞT c dt

ð2Þ

where Mcp and Mb denote the masses of the ceramic plate and cold plate, respectively, and Ccp and Cb the heat capacities; ap and an are the Seebeck coefficients of the P-type and N-type elements; Q is the cooling load absorbed from the external object; Qcp and Qcn represent the heat conduction rates between the ceramic plate and the P-type and N-type elements which can be calculated by

P @T p ðx; y; z; tÞ  dA @x x¼0 N Z @T n ðx; y; z; tÞ ¼ kn  dA @x x¼0

Q cp ¼ Q cn

Z

kp

ð3aÞ ð3bÞ

where kp and kn represent the heat conductivities of the P-type and N-type elements, respectively, and A represents the cross-section area of P-type or N-type elements. In Eqs. (3a) and (3b), it is noted that the heat transfer rate is determined in terms of the integration of the heat flux in direction normal to the plane. 2.2. Hot end (region 2) As also indicated in Fig. 1(a), the hot end region includes a ceramic plate and a heat sink that is used to dissipate heat to the ambient air. The energy equation for the hot end region is written as

ðM s C s þ Mhp C hp Þ

dT h ¼ Iðap  an ÞT h  Q hp  Q hn  hAs ðT h  T 1 Þ dt ð4Þ

where Ms, Mhp, Cs, and Chp are the masses and the heat capacities of the heat sink and the ceramic plate; Qhp and Qhn represent the heat conduction rates between the ceramic plate of hot end and the Ptype and N-type elements which can be calculated by

p @T p ðx; y; z; tÞ  dA @x x¼L N Z @T n ðx; y; z; tÞ ¼ kn  dA @x x¼L

Q hp ¼ Q hn

Z

kp

ð5aÞ ð5bÞ

2004

C.-H. Cheng et al. / International Journal of Heat and Mass Transfer 53 (2010) 2001–2011

2.3. P-type and N-type thermoelectric elements (regions 3 and 4) An energy balance analysis for the P-type element leads to the energy equation: 2

C p qp

2

2

@T p @ Tp @ Tp @ Tp þ þ ¼ kp @t @x2 @y2 @z2

! 

lp @T p A

I

@x

þ

ep A2

I2

ð6aÞ

where Cp, qp, lp, and ep represents the specific heat, the density, the Thomson coefficient, and the electrical resistivity of P-type element, respectively, and I is the electrical current. Similarly, one may have the energy equation for the N-type element as

C n qn

@T n @2T n @2T n @2T n þ þ ¼ kn @t @x2 @y2 @z2

! þ

ln @T n A

I

@x

þ

en A2

I2

ð6bÞ

2.4. Boundary conditions The boundary conditions on the surfaces of the P-type and the N-type elements associated with Eqs. (6a) and (6b) are listed as follows: (1) Surfaces at yp = 0 and yn = 0

h i @T p ðx; y; z; tÞ ¼ ha T p ðx; y; z; tÞjPy¼0  T 1 ; and @y h i @T n ðx; y; z; tÞ kn ¼ ha T n ðx; y; z; tÞjNy¼0  T 1 @y kp

ð7aÞ ð7bÞ

Table 1 Fixed parameters in computation. Parameter

ap an qp qn lp ln ep en

h i @T p ðx; y; z; tÞ ¼ ha T p ðx; y; z; tÞjPz¼0  T 1 ; and @z h i @T n ðx; y; z; tÞ kn ¼ ha T n ðx; y; z; tÞjNz¼0  T 1 @z

kp

where

ð8bÞ

QL ¼ Q

ð9bÞ

(4) Surfaces at zp = Lbp and zn = Lbn

h i @T p ðx; y; z; tÞ ¼ ha T p ðx; y; z; tÞjPz¼Lbp  T 1 ; and @z h i @T n ðx; y; z; tÞ  kn ¼ ha T n ðx; y; z; tÞjNz¼Lbn  T 1 @z

 kp

ð10bÞ

(5) At xp = 0 and xn = 0

T p ðx; y; z; tÞjx¼0 ¼ T c ðtÞ; and

ð11aÞ

T n ðx; y; z; tÞjx¼0 ¼ T c ðtÞ

ð11bÞ

(6) At xp = L and xn = L

T p ðx; y; z; tÞjx¼L ¼ T h ðtÞ; and

ð12aÞ

T n ðx; y; z; tÞjx¼L ¼ T h ðtÞ

ð12bÞ

J kg K J kg1 K1 J kg1 K1 J kg1 K1 J kg1 K1 J kg1 K1 m m kg kg kg kg W m1 K1 W m1 K1 K K K V K1 V K1 kg m3 kg m3 V K1 V K1 Xm Xm

[8] [8] [8] [8] [19–23] [19–23] Measured Measured Measured No cold plate Measured Measured [18] [18] Measured Measured Measured [18] [18] Measured Measured [18] [18] [15] [15]

QL W

ð13Þ

ð14aÞ

  ep 2 en W ¼ Iðap  an ÞðT h  T c Þ þ I V p þ 2 I2 V n 2 A A   lp @T p ln @T n Vp þ I Vn þ  I A @x A @x

ð14bÞ

Note that the net work W is influenced by the Seebeck effect, the Joule heating, and the Thomson effect. The cooling load per unit area (q), in W/cm2, is calculated with

q¼

ð10aÞ

Source 1

The coefficient of performance (COP) is probably the most important index regarding the performance of the thermoelectric coolers. In this study, it is defined by

ð8aÞ

ð9aÞ

1

2.5. Coefficient of performance

COP ¼

(3) Surfaces at zp = 0 and zn = 0

Unit

419 400 900 419 200 200 1.5  103 1.5  103 2.95  105 0 8.64  105 2.95  105 1.8 2.2 300 298 298 0.00015 0.00019 10922.08 10922.08 0.00027 0.000156 1.2  105 1  105

Ccp Cb Cs Chp Cp Cn La Lb Mcp Mb Ms Mhp kp kn Tp T1 Tn

(2) Surfaces at yp = Lap and yn = Lan

h i @T p ðx; y; z; tÞ  kp ¼ ha T p ðx; y; z; tÞjPy¼Lap  T 1 ; and @y h i @T n ðx; y; z; tÞ  kn ¼ ha T n ðx; y; z; tÞjNy¼Lan  T 1 @y

Value

Q Ah

ð15Þ

where Q is the cooling load and Ah is the hot end area, which is fixed at 0.125 cm2 herein. The magnitude of the coefficient of performance (COP) of the thermoelectric cooler is calculated when the steady regime is attained. The transient variation in temperatures and the COP value at steady state are predicted in wide ranges of physical and geometrical parameters. The varied parameters in this study include

Table 2 Specific heats of thermoelectric materials. Material

C (250 K)

C (298 K)

C (400 K)

Unit

Reference

Bi2Te3 Bi2Te3 Bi2Se3 Sb2Te3 Bi0.52 Sb1.48Te3

186.587 154.412 197.218 203.257 182.860

189.834 157.588 203.613 205.493 187.059

194.246 163.285 208.774 209.963 194.473

J kg1 K1 J kg1 K1 J kg1 K1 J kg1 K1 J kg1 K1

[19] [20] [21] [22] [23]

C.-H. Cheng et al. / International Journal of Heat and Mass Transfer 53 (2010) 2001–2011

with all the cases in computation. Densities of the P-type and Ntype elements (10922.08 kg m3) are determined by direct measurement in terms of the mass and volume of the element. And, according to the measurement data, no appreciable difference in density between these two elements is found. Note that the values of the specific heats, Cp and Cn, are dependent on the thermoelectric material used; however, it is found that even for the same material, the information is rather scattering. Some of the existing information of the specific heat taken from Refs. [19–23] are listed in Table 2. It can be found that the specific heat of the major thermoelectric materials is within 154–210 J kg1 K1 and according to the existing data, no appreciable difference between the P-type and N-type materials is found. Therefore, in this study, the values of the specific heats are assigned to be Cp = 200 J kg1 K1 and Cn = 200 J kg1 K1.

Mesh 330

11x11x16 21x21x31 31x31x46 41x41x61

Th

315

T/(K)

2005

Tc 300

3. Numerical methods

285 100

200

300

400

500

600

700

Time/(s) Fig. 2. Grid-independency check by a comparison in hot end and cold end temperatures obtained from different grids.

Table 3 Base-line case parameters. Parameter

Value

Unit

Ah As h ha I L q

1.25  105 4.254  105 20 10 1 2.325  103 800

m2 m2 W m2 K1 W m2 K1 A m W m2

the applied current (I), total heat transfer surface area of the heat sink (As), heat transfer coefficient on surface of the heat sink (h), heat transfer coefficient on side surfaces of thermoelectric elements (ha), height of the thermoelectric elements (L), and cooling load per unit area (q). On the other hand, the parameters are determined either according to direct measurement or from the existing information [8,15,18–23]. Table 1 displayed the fixed parameters

The governing equations are discretized based on the finite-difference method. These time derivative terms in Eqs. (6a) and (6b) are treated by the forward-difference scheme. A 31  31  46 grid point system is adopted for the solution domain y  z  x after a careful grid-independence check. Fig. 2 displays the grid-independence of the solutions for the case with As = 4.254  105 m2, h = 20 W m2 K1, ha = 10 W m2 K1, I = 1 A, L = 2.25  103 m, and q = 800 W cm2. Other fixed parameters are equal to those given in Table 3. In this figure, four grid systems, 11  11  16, 21  21  31, 31  31  46, and 41  41  61, are tested, and the comparison in the time histories of the hot end and the cold end temperatures for the case are presented. In addition, the magnitudes of COP by different grids are also provided in a small table shown in the figure. It is found that for this case the relative errors in the temperature solutions and COP between 31  31  46 and 41  41  61 grids are all within 0.1%. However, note that the computation time required with the 41  41  61 grids is three times that with the 31  31  46 grids. Hence, to save computation efforts without loss in accuracy, the 31  31  46 grids are adopted. 4. Experiments By conducting experiments for the thermal behavior of the thermoelectric module, the validity of the simulation model is possible to be confirmed in part. Herein a zero-cooling-loading test has

Heat sink Thermoelectric cooler module Ceramic fiber blanket Thermocouples

Data logger Power supply Fig. 3. Experimental system.

PC

2006

C.-H. Cheng et al. / International Journal of Heat and Mass Transfer 53 (2010) 2001–2011

been performed for validation of the simulation model. The reason for performing the zero-cooling-load test is to avoid the errors in measuring the magnitude of the cooling load and the heat losses with it which could lead to uncertainties in the testing conditions. The experimental apparatus of the experiments is shown in Fig. 3, which consists of the following components:

(2) Heat sink: A plate-fin heat sink is attached on the top of the ceramic plate at hot end for heat dissipation to the ambient. The heat sink is made of aluminum 6063 with a total heat transfer area of As = 4.254  105 m2. Total heat transfer area of the heat sink was divided by the number of element pairs to yield the heat transfer area per element pair, which is used in the one pair simulation. (3) Ceramic fiber blanket: To fulfill a zero-cooling-loading condition, a ceramic fiber blanket (Al2O3 and SiO2) is placed under the cold end of the thermoelectric model for insulation. (4) Power supply: The power supply is used to supply a DC current pass through the device. The accuracy of the power supply is ±0.5%. (5) Thermocouples: Three T-type thermocouples (copper/copper-constantan alloy) are installed for temperature

(1) Thermoelectric cooler module: The thermoelectric cooler module used is a typical one which comprises 127 pairs of P-type and N-type elements connected in series which are clamped by two ceramic plates. The area of the thermoelectric cooler module is 4 cm4 cm. For this particular case, the dimensions of the P-type and the N-type elements are L = 2.325  103 m, La = 1.5  103 m, and Lb = 1.5  103 m. In the zero-cooling-load case, q = 0.

320

360

I = 0.5A

Th

I = 1A

Th

350 340

T/(K)

310

T/(K)

330

Tc

300

Tc

320 310 300

-Numerical - Experimental

290

- Numerical - Experimental

290 280

0

100

200

300

400

0

500

100

200

300

Time/(s)

Time/(s)

(a) I = 0.5 A

(b) I = 1 A

400

500

425

I = 1.5A

Th

400

375

T/(K)

Tc 350

325

300

- Numerical - Experimental

275 0

100

200

300

400

500

Time/(s)

(c) I = 1.5 A Fig. 4. Comparison in temperature variation between numerical predictions and experiments at zero-cooling-load per unit area (q = 0). (a) I = 0.5 A, (b) I = 1 A, and (c) I = 1.5 A.

2007

C.-H. Cheng et al. / International Journal of Heat and Mass Transfer 53 (2010) 2001–2011

measurement, whose error is within ±0.5°. For monitoring the hot and cold ends temperatures, one thermocouple is placed at the interface between the hot end and the heat sink, and one at the interface between the cold end and the ceramic fiber blanket. The third thermocouple is installed for measuring the ambient temperature. (6) Data logger: Fluke 2620A Data Logger is used for temperature data collection. An RS232 serial interface is used to connect the data logger to a personal computer for real time data acquisition. The experiment is started as the DC current is applied steadily by the power supply. The temperature data of the hot and the cold sides measured by the thermal couples are transmitted to data logger and transmitted into computer via RS232 for analysis. The

I [A] 2 1.5 1 0.5

360

Th

Th Th

Th/(K) & Tc/(K)

As [m2]2 As [m ] 8.508×10-5 8.508E-5 4.254×10-5 4.254E-5-5 2.127×10 2.127E-5

330

Th

340

Tc

320

Th Tc

Th 300

Th/(K) & Tc/(K)

380

geometrical parameters of the simulation model are identical to the specifications of the experimental device. The comparisons between the numerical predictions and experiments are made at I = 0.5, 1.0, and 1.5 A. All the fixed parameters in the computation are already given in Table 1. Note that the heat transfer coefficients h and ha are assigned to be 20 and 10 W m2 K1, respectively. The values are corresponding to the air natural convection situation based on the operating conditions of experiments. Fig. 4 shows the comparison in temperature variation between the numerical predictions and the experiments at q = 0. It is found by experiments that the hot end temperature increases with time monotonically. For the case at I = 0.5 A, the hot end temperature reaches a temperature of 313 K in 500 s; however, in the same period of time it reaches 422 K for the case at I = 1.5 A. On the other hand, the cold end temperature descends

315

Th Tc 300

Tc

Tc

Tc

Tc 285

280

100

200

300

400

500

0

100

200

300

400

500

Time/(s)

Time/(s)

(a) Th and Tc

Fig. 6. Effects of total area of the heat sink on transient temperature variation for base-line case at steady state.

40

h

340

Th

[W/K-m2]

40 20 5 Th Tc Th

20

Th/(K) & Tc/(K)

COP

320

300

Tc Tc 280

0 0

0.5

1

1.5

I/(A)

(b) COP Fig. 5. Effects of applied current on transient temperature variation and COP for base-line case at steady state. (a) Th and Tc; (b) COP.

0

100

200

300

400

500

Time/(s) Fig. 7. Effects of heat transfer coefficient on the surface of the heat sink on transient temperature variation for base-line case at steady state.

2008

C.-H. Cheng et al. / International Journal of Heat and Mass Transfer 53 (2010) 2001–2011

from start to a minimum value and then goes up to a value higher than the ambient temperature. This may be explained by Eq. (2). At beginning, the Seebeck effect causes heat output from the cold end so that the cold end temperature is decreased. However, after a period of time, when the temperature of hot end is increased, the heat conduction from the hot to the cold end is also increased and it results in the increase in the cold end temperature. The difference between the hot- and the cold end temperatures is approximately maintained at 12, 25, and 40 K for the cases at I = 0.5, 1.0, and 1.5 A, respectively. The numerical predictions by the simulation model closely agree with the experimental data. Hence, the validity of the present simulation model can be ensured.

parameters. In this study, for comparisons a base-line case is chosen and analyzed firstly. The physical and geometrical parameters with the base-line case are displayed in Table 3. It is noted that the parameters of the base-line case are identical to the conditions of the experiments described in the subsequent section. Also note that the current with the base-line case is assigned to be 1 A according to the specifications of the existing commercial TEC modules. In the parametric study, when a parameter is varied to evaluate its effects, other parameters remain the same as with the base-line case. The following figures are plotted for the baseline case; therefore, the conditions for the case are already given in Tables 1–3.

5. Results and discussion

5.1. Effects of applied current

Next, the simulation model is adopted to perform a parametric study to investigate the effects of the physical and geometrical

Fig. 5 shows the effects of applied current on the transient variation in temperatures and the COP value at steady state for the 380

15 10 5

Th/(K) & Tc/(K)

320

L

[W/K-m2]

Th

300

Tc

[mm]

10 5 2.325

360

Th/(K) & Tc/(K)

ha

Th

340

Th

320

Th

300

Tc Tc Tc

280

280 260

0

100

200

300

400

0

500

200

300

Time/(s)

(a) Th and Tc

(a) Th and Tc

0.49

400

500

0.5

0.489

0.45

0.488

0.4

0.487

0.35

COP

COP

100

Time/(s)

0.486

0.3 0.25

0.485

0.2 0.484

0.15 0.483

0.1 10

20 /(Wm-2K-1)

30

ha

(b) COP Fig. 8. Effects of heat transfer coefficient on the surfaces of the P-type and N-type elements on transient temperature variation and COP for base-line case at steady state. (a) Th and Tc; (b) COP.

2

3

4

5

6

7

8

9

L/(mm)

(b) COP Fig. 9. Effects of height of the P-type and N-type elements height of the P-type and N-type elements on transient temperature variation and COP for base-line case at steady state. (a) Th and Tc; (b) COP.

2009

C.-H. Cheng et al. / International Journal of Heat and Mass Transfer 53 (2010) 2001–2011

cases with I = 0.5, 1.0. 1.5, and 2.0 A. Note that for this case, the cooling load per unit area q is fixed at 800 W cm2. As shown in this figure, the temperature of the hot end and the temperature difference between the hot and cold ends both increase with the applied current. This is because in Eq. (6a), the heat sink term due to Thomson effect is proportional to current while the heat source term due to Joule heating is proportional to the square of current. As a result, once the value of current is elevated, the heat source term gradually dominates, and hence the slope of increase in the temperatures of the P-type and the hot end are both higher. In a similar manner, even the cold end temperature is decreased at the beginning, it is increased to a value exceeding the ambient temperature in the steady state regime. However, if the magnitude of the current is not high enough (say, 0.5 or 1.0 A), the cold end temperature will still be maintained at a lower value. The COP is displayed in Fig. 5(b) as a function of applied current. It is found that the magnitude of COP is rapidly decreased as the current is increased. This may be attributed to the increase in the net work (W), which is defined by Eq. (14b), resulting from an increase in current.

0.01 m, and it is found that the temperature difference between the two ends increases with L. However, as L is changed from 0.002325 to 0.01 m, the value of COP is reduced approximately from 0.5 to 0.1. 5.6. Effects of cooling load The cooling load could be different in response to various situations. Under different cooling loads, the thermal behavior of the thermoelectric cooler could be different. Effects of cooling load on thermal behavior are illustrated in Fig. 10 for the base-line case. In the cases shown in Fig. 10(a), the value of q is assigned to be 0, 400, 800, and 2400 W cm2. It is observed that both the temperatures of the hot and the cold ends increase with q. The temperature difference (Th  Tc) reaches a maximum value of 24.51 K for the case at q = 0. Furthermore, as shown in Fig. 10(b), the value of COP linearly increases with q. This complies with Eq. (13) in which the COP appears to be proportional to the cooling load per unit area 2

q [W/cm ] 0.24 0.08 0.04 0

5.2. Effects of total area of the heat sink

340

Th

Th/(K) & Tc/(K)

Convective heat transfer from the heat sink to the ambient could be influential. Fig. 6 shows the results for the cases with As varied from 2.127  105 to 8.508  105 m2. According to Eq. (4) for Th, it is expected that the increasing rate of hot end temperature, dTh/dt, is lowered when the magnitude of hAs(Th  T1) is increased. The curves plotted in Fig. 6 agree with the expectation. The hot end temperature is remarkably reduced by increasing the total area of the heat sink. Note that the dependence of cold end temperature on total area of the heat sink is similar to that of the hot end temperature due to the conduction through the semiconductor elements between the two ends.

Th

Th

Th

320

Tc

Tc

300

Tc

Tc

5.3. Effects of heat transfer coefficient on the surface of the heat sink 280

Fig. 7 shows the effects of heat transfer coefficient on the surface of the heat sink on transient temperature variation. Basically, to increase heat dissipation by convection one may increases the value of either As or h, if the temperature difference Th  T1 is kept constant. The effect of h on the thermal behavior of thermoelectric cooler is similar to that of As. When the value of h is assigned to be 5, 20, and 40 W m2 K1, the steady state hot end temperature is 340.9 K, 318.8 K, and 309.9 K and the cold end temperature 313.7 K, 294.3 K, and 286.4 K, respectively.

100

200

300

400

500

Time/(s)

(a) Th and Tc 5 4.5 4

5.4. Effects of heat transfer coefficient on the surfaces of the P-type and N-type elements

5.5. Effects of height of the P-type and N-type elements Effects of the height of the semiconductor elements are exhibited in Fig. 9. The influence of an increase in element height is actually involved because it increases the thermal resistance against heat conduction from the hot end to the cold but also increases electrical resistance as well as the Joule heating effects. In Fig. 9(a), the value of L is assigned to be 0.002325, 0.005 m, and

COP

A heat loss from the surfaces of the P-type and N-type elements causes a reduction in the COP. Fig. 8 shows the effects of ha on transient temperature variation and COP at steady state. In Fig. 8(b), it is observed that as ha is elevated from 5 to 15 W m2 K1, the value of COP is reduced from 0.4896 to 0.483. Meanwhile, an increase in ha results in a decrease in both Th and Tc; therefore, the temperature difference is only slightly influenced by a change in the heat transfer coefficient.

3.5 3 2.5 2 1.5 1 0.5 0

0

0.2

0.4

0.6

0.8

q/(Wcm-2)

(b) COP Fig. 10. Effects of cooling load on transient temperature variation and COP for base-line case at steady state. (a) Th and Tc; (b) COP.

2010

C.-H. Cheng et al. / International Journal of Heat and Mass Transfer 53 (2010) 2001–2011

q as long as the magnitude of the net work W is fixed. In this case, the temperature difference between the hot and the cold ends is only slightly decreased by increasing the value of q. Therefore, the magnitude of the net work is not appreciably altered. In addition, the Joule heating and the Thomson terms in Eqs. (6a) and (6b) can be regarded as a heat source and a heat sink terms, respectively. Therefore, in general the Joule heating leads to a reduction in the magnitude of COP, whereas the Thomson effect tends to increase the COP. However, the investigation of the subtle influence of these effects is still open to discussion and worth further investigation. As stated earlier, the present three-dimensional model treats the P-type and the N-type semiconductor elements as two separate parts. Therefore, one may take advantages of the three-dimen-

sional model particularly when the properties of the two elements are rather different. Plotted in Fig. 11 are the planar temperature distributions at different z- and x-locations within the P-type and the N-type elements, of which the electrical resistivities are set to be 3  105 and 6.5  106 X m, respectively. Note that these two values are still within the reasonable range of electrical resistivity in practical cases [15]. Other parameters are equal to those of the base-line case. Here, the difference in electrical resistivity between the two elements results in the discrepancy in temperature distribution. The steady state temperature contours on the x–y cross-sectional planes at z = 0, Lb/4, and Lb/2, and on the y–z cross-sectional planes at x = L/4, L/2, and 3L/4, are plotted and the discrepancy in temperature distributions between the two elements is clearly seen.

Fig. 11. Steady-state temperature distributions at different locations within the P-type and the N-type elements. (a) At different z-locations. (b) At different x-locations.

C.-H. Cheng et al. / International Journal of Heat and Mass Transfer 53 (2010) 2001–2011

6. Conclusions In this model, the solution domain is divided into four major regions, namely, cold end (region 1), hot end (region 2), and the Ptype and N-type thermoelectric elements (regions 3 and 4). The present model treats the P-type and the N-type semiconductor elements as two separate parts so that the difference in thermal behavior between the two elements can be evaluated. To verify the numerical predictions, experiments have been conducted to measure the temperature variations in the cold and the hot ends for a certain particular cases. Close agreement between the numerical and the experimental data has been observed. The simulation model is adopted to perform a parametric study to investigate the effects of the parameters such as electrical current (I), cooling load per unit area (q), heat transfer coefficients on the surfaces (h and ha) and geometrical parameters such as total area of heat sink (As) and height of semiconductor elements (L). The following conclusions are yielded based on the obtained results: (1) The temperature of the hot end and the temperature difference between the hot and cold ends both increase with the applied current. It is also found that the magnitude of COP is rapidly decreased as the current is increased. This may be attributed to the increase in the net work input resulting from an increase in current. (2) Convective heat transfer from the heat sink to the ambient could be influential. The hot end temperature is remarkably reduced by increasing the convective heat transfer of the heat sink. The dependence of cold end temperature on the heat rejection from heat sink is similar due to the conduction through the semiconductor elements between the two ends. (3) A heat loss from the surfaces of the P-type and N-type elements causes a reduction in the COP; however, the temperature difference is only slightly influenced by a change in the heat transfer coefficient. (4) The influence of an increase in element height is actually involved. It is found that the temperature difference between the two ends increases with L. However, as L is changed from 0.002325 to 0.01 m, the value of COP is reduced approximately from 0.5 to 0.1. (5) It is observed that both the temperatures of the hot and the cold ends increase with cooling load q, and the value of COP linearly increases with q. In the cases considered, the temperature difference between the hot and the cold ends is only slightly decreased by increasing the value of q. Acknowledgements The authors would like to thank the National Science Council, Taiwan, ROC, for their financial support under Grant: NSC 97-

2011

2221-E-006-111-MY3. Authors would also like to thank the National Space Organization, Taiwan, ROC, for their partial financial support and help in experiments.

References [1] Y. Kraftmakher, Simple experiments with a thermoelectric module, Euro. J. Phys. 26 (2005) 959–967. [2] G.D. Mahan, Figure of merit for thermoelectrics, J. Appl. Phys. 65 (1989) 1578– 1583. [3] R. Venkatasubramanian, E. Siivola, T. Colpitts, B. O’Quinn, Thin-film thermoelectric devices with high room-temperature figures of merit, Nature 413 (2001) 597–602. [4] L.S. Stilbans, A.R. Regel, Thermoelect. Power Generat. 1 (1968) 1341. Soviet Physics Semiconductors-USSR. [5] J.E. Parrott, A new theory of size effect in electrical conduction, Proc. Phys. Soc. Lond. 85 (1965) 1143. [6] P.E. Gray, The Dynamic Behavior of Thermoelectric Devices, Wiley, New York, 1960. [7] R.P. Bywaters, H.A. Blum, The transient behavior of cascade thermoelectric heat pumps, Energy Convers. 10 (1970) 193–200. [8] B.J. Huang, C.L. Duang, System dynamic model and temperature control of a thermoelectric cooler, Int. J. Refrigerat. 23 (2000) 197–207. [9] B.J. Huang, C.L. Duang, A design of thermoelectric cooler, Int. J. Refrigerat. 23 (2000) 208–218. [10] M. Naji, M. Alata, M.A. Al-Nimr, Transient behavior of a thermoelectric device, Proc. Inst. Mechan. Eng. A: J. Power Energy 217 (2003) 615– 621. [11] X.C. Xuan, K.C. Ng, C. Yap, H.T. Chua, A general model for studying effects of interface layers on thermoelectric devices performance, Int. J. Heat Mass Transfer 45 (2002) 5159–5170. [12] M.J. Huang, R.H. Yen, A.B. Wang, The influence of the Thomson effect on the performance of a thermoelectric cooler, Int. J. Heat Mass Transfer 48 (2005) 413–418. [13] R.A. Khire, A. Messac, S.V. Dessel, Design of thermoelectric heat pump unit for active building envelope systems, Int. J. Heat Mass Transfer 48 (2005) 4028– 4040. [14] A. Chakraborty, B.B. Saha, S. Koyama, K.C. Ng, Thermodynamic modeling of a solid state thermoelectric cooling device: temperature-entropy analysis, Int. J. Heat Mass Transfer 49 (2006) 3547–3554. [15] K.H. Lee, O.J. Kim, Analysis on the cooling performance of the thermoelectric micro-cooler, Int. J. Heat Mass Transfer 50 (2007) 1982–1992. [16] A.J. Minnich, M.S. Dresselhaus, Z.F. Ren, G. Chen, Bulk nanostructured thermoelectric materials: current research and future prospects, Energy Environ. Sci. 2 (2009) 466–479. [17] W.E. Beadle, J.C. Tsai, R.D. Plummer, Quick Reference Manual for Silicon Integrated Circuit Technology, John Wiley& Sons, 1985. [18] D.M. Rowe (Ed.), CRC Handbook of Thermoelectric, CRC Press, USA, 1995, pp. 145–230. [19] D.M. Rowe (Ed.) Thermoelectrics Handbook, CRC Press, USA, 2006, pp. 27–28. [20] N.P. Gorbachuk, A.S. Bolgar, V.R. Sidorko, L.V. Goncharuk, Heat capacity and enthalpy of Bi2 Se3 and Bi2 Te3 in the temperature range 58–1012 K, Powder Metallur. Metal Ceram. 43 (2004) 284–290. [21] A.S. Pashinkin, A.S. Malkova, Heat capacity of solid bismuth selenide (Bi2 Se3), Zhurnal Fizicheskoi Khimii 79 (2005) 1325–1327. [22] A.S. Pashinkin, A.S. Malkova, M.S. Mikhailova, Heat capacity of solid antimony telluride, Zhurnal Fizicheskoi Khimii 81 (2008) 1–3. [23] Y.I. Shtern, A.S. Malkova, A.S. Pashinkin, V.A. Fedorov, Heat Capacity of the nBi2Te2.88Se0.12 and p-Bi0.52Sb1.48Te3 Solid Solutions, Inorgan. Mater. 44 (2008) 1057–1059.