Applied Energy 108 (2013) 340–348

Contents lists available at SciVerse ScienceDirect

Applied Energy
journal homepage: www.elsevier.com/locate/apenergy

Transient modeling and dynamic characteristics of thermoelectric cooler

Jing-Hui Meng a, Xiao-Dong Wang b,c,⇑, Xin-Xin Zhang a
a
School of Mechanical Engineering, University of Science and Technology Beijing, Beijing 100083, China
b
State Key Laboratory of Alternate Electrical Power System with Renewable Energy Sources, North China Electric Power University, Beijing 102206, China
c
Beijing Key Laboratory of Multiphase Flow and Heat Transfer for Low Grade Energy, North China Electric Power University, Beijing 102206, China

h i g h l i g h t s

 A complete three-dimensional TEC transient model was proposed.

 The model couples heat conduction and electric conduction.
 Dynamic behaviors are studied at extensive operating conditions.
 Temperature-dependent material properties have strong effect on dynamic behaviors.

a r t i c l e i n f o a b s t r a c t

Article history: Dynamic characteristics are extremely important for design and operation of thermoelectric coolers
Received 18 October 2012 (TECs). This paper develops a three-dimensional transient TEC model based on the coupling of heat trans-
Received in revised form 10 January 2013 fer and electric conduction within semiconductors. The model takes into account all thermoelectric
Accepted 22 March 2013
effects, including Joule heating, Thomson effect, Peltier effect and Fourier’s heat conduction. For most
Available online 11 April 2013
of semiconductor materials, Seebeck coefficient, electric conductivity and thermal conductivity are
strongly temperature-dependent. Therefore, the present transient model is used to compare dynamic
Keywords:
temperature variations at the cold and hot ends with constant and variable material properties. Small,
Thermoelectric cooler
Dynamic characteristics
medium, and large applied currents with various cooling loads are adopted as operating conditions.
Heat transfer The results show that, at small currents, constant property model developed by this work can predict
Peltier accurately the dynamic characteristics, however, with the increase in current, the temperature-depen-
Seebeck dence of properties have more and more remarkable effect on the dynamic temperature variations, espe-
Thomson cially for high cooling loads. When the current is larger than a specific value, the heat transferred from the
hot end to the cold end by Fourier’s heat conduction will exceed the heat adsorbed at the cold end by
Peltier effect, thus, the temperatures at the cold and hot ends increase continuously, the TEC cannot reach
the steady-state. This phenomena is predicted exactly by the variable property model, oppositely, the
constant property model predicts that the TEC still supply refrigeration.
Ó 2013 Elsevier Ltd. All rights reserved.

1. Introduction ments are composed of a pair of p- and n-type semiconductors.

When electric current flows across the thermoelectric elements,
In recent years, TECs have attracted more and more attention the heat is transported from the cold end to the hot end caused
due to their features of no compressors and refrigerants, large by the Peltier effect.
operating temperature range, easy to control, reliable operation, A large amount of works considered the TEC performance
layout flexibility, adaptability and other characteristics [1–3]. The [5–13] and the TEG performance [14–22] under steady-state
TEC has been frequently used for the cooling of electric devices condition, while few works [4,23–29] investigated its transient or
such as CPU, infrared sensor, ice-point reference in thermocouple dynamic characteristics. Actually, transient characteristics are also
thermometry, and refrigerators [4]. Generally, the TEC includes a extremely important in practical TEC operation. First, cooling load
number of thermoelectric elements, which are connected electri- at the cold end and/or cooling capacity of heat sink and ambient
cally in series and thermally in parallel. The thermoelectric ele- temperature at the hot end may be variable [24]. Second, the start
up and shut down characteristics of TECs are also a major concern
for cooling purpose [25,26]. For example, for a fixed cooling load at
⇑ Corresponding author at: State Key Laboratory of Alternate Electrical Power
the cold end and a fixed cooling capacity of heat sink at the hot
System with Renewable Energy Sources, North China Electric Power University,
Beijing 102206, China. Tel./fax: +86 10 62321277. end, when a specific current is applied to the TEC, how much is
E-mail address: [email protected] (X.-D. Wang). the consuming for the TEC to reach the steady-state and what

0306-2619/$ - see front matter Ó 2013 Elsevier Ltd. All rights reserved.
http://dx.doi.org/10.1016/j.apenergy.2013.03.051
 J.-H. Meng et al. / Applied Energy 108 (2013) 340–348 341

are the corresponding steady-state temperatures of cold and hot The heat is adsorbed at the cold end and liberated at the hot end
ends. In addition, as a current pulse with a magnitude of several caused by the Peltier effect, forming a temperature difference
times higher than the steady-state optimum current is applied to DT = TH  TL. With the temperature difference, the heat will be
a TEC, the transient lower temperature than that reachable at the transferred from the hot end to the cold end due to Fourier’s heat
steady-state can be achieved at the cold end due to the delay of conduction. In the role of the current and temperature gradient,
the thermal diffusion of the volumetric Joule heat, which is re- the Joule heat and the Thomson heat will be generated within
ferred to as the transient thermoelectric effect [27,28]. semiconductors. Thus the final cooling capacity of the TEC, QL, is
Up to now, some transient models have been developed to determined by cooperative effect of Peltier heat, Fourier’s heat
understand the dynamic behaviors of TEC. Mostly the previous tran- conduction, Joule heat, and Thomson heat.
sient models were limited to one-dimensional problems [4,24,28],
where the Thomson effect was ignored and constant material prop- 2.1. Governing equations
erties were used because small temperature difference between the
cold and hot ends was assumed. In addition, as mentioned by Cheng The transient TEC model is developed from our previous steady-
et al. [26], in the existing models, the p-n element pair was simply state model [29] by adding a transient term. The basic governing
treated as a single bulk material so that the difference in thermal equations include the energy equations and the electric potential
behavior between the two semiconductor elements was not possi- equations. The energy equations of connectors, p- and n-type semi-
ble to evaluate. Hence, Cheng et al. developed a three-dimensional conductors, ceramic plates, and heat sink are as follows:
model with p-type and n-type semiconductors as two separate
parts and the Thomson effect was taken into account. However,
@T J2
the material properties were still assumed to be constants in the ðqcp Þconn ¼ r  ðkconn rTÞ þ b ~ J  rT ð1Þ
Cheng’s model.
@s rconn conn
Recently, we proposed a general, three-dimensional TEC steady-
state model [29], which is different from previous models [5–13] in
which only heat conduction equation with Joule heat and/or Thom- @T J2
son heat as internal heat sources is solved. Our model introduced ðqcp Þp ¼ r  ðkp rTÞ þ b~
J  rT ð2Þ
@s rp p
the coupling of heat conduction and electrical conduction, and con-
sidered all the effects occurred in the TEC. The model was used to
@T J2
figure out the performance of TECs with the temperature-depen- ðqcp Þn ¼ r  ðkn rTÞ þ b~
J  rT ð3Þ
dent material properties. The predictions showed that the variable @s rn n
properties and the heat losses to the ambient gas have significant
effects on the cooling capacity and the coefficient of performance
(COP) of the TEC. Three-dimensional temperature distributions @T
within semiconductors was observed and it became more remark- ðqcp Þcer ¼ r  ðkcer rTÞ ð4Þ
@s
able at large temperature differences and high currents.
The purpose of this work reaches the following targets: (1) @T
extending our previous steady-state model [29] to a transient ðqcp Þsink ¼ r  ðksink rTÞ ð5Þ
@s
one; (2) employing the developed transient model to investigate
the dynamic characteristics of TECs under various operating condi- where q is the density, cp is the specific heat, k is the thermal
tions. The dynamic behaviors for TECs with temperature-dependent conductivity, r is the electric conductivity, and b is the Thomson
properties are analyzed and compared with those constant coefficient. The subscripts conn, p, n, cer, and sink denote connector,
properties. p- and n-type semiconductors, ceramic plate, and heat sink, respec-
tively. The first terms on the right side in Eqs. (1)–(5) denote Fou-
rier’s heat conduction, the second and third terms denote internal
2. Model development heat sources due to Joule heating and Thomson effect, respectively.
J(x, y, z) is the local current density, which is assumed to be constant
The schematic of a TEC device containing one semiconductor and equal to I/A in the previous models [4–13,23–28], where I is the
element pair (referred to as TEC element hereafter) is shown in total applied current, and A is the cross-sectional area of semicon-
Fig. 1. A TEC device consists of a number of TEC elements. In con- ductors. The Thomson coefficient can be derived from Seebeck coef-
sideration that the thermal characteristics among the TEC ele- ficient, or:
ments are periodic, hence, only one element is considered here.
The TEC element includes a pair of p- and n-type semiconductor da
b¼T ð6Þ
columns, three metallic connectors, and two electrically insulating dT
ceramic plates. The ceramic plate, connector, and semiconductor where a is the Seebeck coefficient. The electric potential is the driv-
columns have thicknesses of H0, H1, and H2, respectively. The semi- ing force of the electrons and holes in the semiconductor, which can
conductor columns have square cross-section with the side length be obtained by solving the following equation:
of L2. The distance between p- and n-type semiconductors is L1. A
heat sink with the thickness of Hs is attached on the hot end of the r  ðrðr/  arTÞÞ ¼ 0 ð7Þ
TEC element for heat dissipation. The geometric parameters of the where / is the electric potential, arT is Seebeck electromotive force
TEC element simulated in the present work are as follows, coming from Seebeck effect. Once / is determined, the electric field
L1 = 0.2 mm, L2 = 0.5 mm, H0 = 0.2 mm, H1 = 0.1 mm, H2 = 1.0 mm, can be calculated by the following equation:
and Hs = 0.2 mm.
When applied current flows from the p-type semiconductor to ~
E ¼ r/ þ arT ð8Þ
n-type semiconductor, the holes in the p-type semiconductor and
Finally, the current density vector in Eqs. (1)–(3) can be calcu-
the electrons in the n-type semiconductor migrate from the cold
lated as follows:
end to the hot end, the corresponding Peltier heats will be gener-
ated at the interface between connectors and semiconductors.
342 J.-H. Meng et al. / Applied Energy 108 (2013) 340–348

Fig. 1. The schematic of TEC element.

~
J ¼ r~
E ð9Þ On the internal interfaces between adjacent materials, the tem-
perature and heat flux are assumed to be continuous. The electric
boundary conditions for the TEC element are as follows: a fixed
2.2. Initial conditions and boundary conditions current is applied to the TEC element and a zero electric potential
is specified on the side surface of connector with x = 2L1 + 2L2
The initial conditions are as follows. The temperature of the TEC (Fig. 1) or:
element is uniform and equals to the ambient temperature of
T = 300 K, and the initial electric potential of the TEC element is 0. Ijx¼0 ¼ constant 0 < y < L2 ; H0 < z < H0 þ H1 ð13Þ
The adopted boundary conditions are: on the top of the heat
sink, heat is transferred to the ambient by convective heat transfer, /jx¼2L1 þ2L2 ¼ 0 0 < y < L2 ; H0 < z < H0 þ H1 ð14Þ
or
On the other surfaces, the current cannot flow out of the TEC
@T element, thus, we have:
ksink ¼ heff ðTjz¼2ðH0 þH1 ÞþH2 þHs  T 1 Þ ð10Þ
@z z¼2ðH0 þH1 ÞþH2 þHs
~
J ~
n¼0 ð15Þ
where heff is the equivalent convective heat transfer coefficient,
which will be determined in the next section.
On the bottom of the TEC element, a constant heat flux is as- 2.3. Material properties
sumed as follows:
The n-type and p-type semiconductors are assumed to be Bi2(-
qL jz¼0 ¼ constant ð11Þ
Te0.94Se0.06)3 and (Bi0.25Sb0.75)Te3. Minnich et al. [3] proposed that
On the side surfaces of the TEC element, the adiabatic boundary the properties of some thermoelectric materials are strongly tem-
conditions are specified, or: perature-dependent, hence for the purpose of constructing a gen-
eral TEC model, the variable properties should be considered. In
@T
¼0 ð12Þ the present work, the thermal conductivity, electric conductivity
@n
and Seebeck coefficient are all assumed to be the functions of tem-
perature as follows:

Table 1
Parameters in Eqs. (16)–(18) for thermoelectric materials (T0 = 300 K).

Material Thermal conductivity Electric resistivity Seebeck coefficient

k0 (W m1 K1) A1 (K1) A2 (K2) q0 (X m) B1 (K1) B2 (K2) a0 (V K1) C1 (K1) C2 (K2)
3 5 6 3 6 4 3
p-Type 1.472 1.29  10 1.35  10 8.826  10 5.88  10 8.93  10 2.207  10 1.55  10 3.15  106
n-Type 1.643 9.80  104 1.56  105 8.239  106 4.70  103 2.67  106 2.230  104 5.62  104 4.65  106
Copper 400 0 0 1.700  109 0 0 6.500  106 0 0
 J.-H. Meng et al. / Applied Energy 108 (2013) 340–348 343

Table 2 Table 3
Properties for AlN ceramic plate (T = 300 K) [22]. Properties of heat sink system.

Parameter Value Unit Parameter Value Unit Source

3
qcer 3250 kg m Fin
kcer 260 W K1 m1 Height Hf 14.0 mm [4]
cp,cer 7.6  105 J kg1 K1 Thickness tf 0.2 mm [4]
Length Lsink 0.5 mm -
Fin to fin spacing Dg 0.2 mm [4]
Base width Wsink 1.4 mm –
Base height bf 0.2 mm [4]
Coolant (air)
Thermal conductivity ka 26.3  103 W m1 K1 [23]
Density qa 1.1614 kg m3 [23]
Viscosity la 184.6  107 N s m2 [23]
Specific heat cp,a 1007 J kg1 K1 [23]

Q H ¼ heff AðT top

sink  T 1 Þ ð19Þ

where QH is the total heat dissipation of real heat sink, T top

sink is the
top temperature of equivalent plate heat sink, and A is the base area
of heat sink, which equals the area of ceramic plate, or
A = Lsink  Wsink = L2  (2L1  2L2).
For a finned structure heat sink, Mereu et al. [31] measured
experimentally the relationship between the heat dissipation with
Fig. 2. The schematic of heat sink. heat sink structure and operating conditions, expressed as:

kðTÞ ¼ kðT 0 Þ½1 þ A1 ðT  T 0 Þ þ A2 ðT  T 0 Þ2  ð16Þ qa W sink Hf D2g Dp

QH ¼ cp;a ðT max  T 1 Þ ð20Þ
1 þ Df
t
12la Lsink
g
1
¼ qðTÞ ¼ qðT 0 Þ½1 þ B1 ðT  T 0 Þ þ B2 ðT  T 0 Þ2  ð17Þ
rðTÞ where qc, lc, cp,c are the density, viscosity, and specific heat of air,
Dp is the pressure drop across the heat sink, Tmax is the maximum
aðTÞ ¼ aðT 0 Þ½1 þ C 1 ðT  T 0 Þ þ C 2 ðT  T 0 Þ2  ð18Þ temperature observed in the heat sink. Copper is adopted as the
material of heat sink, and its thermal conductivity is 400 W m1 K1
where T0 = 300 K is the reference temperature, k(T0), q(T0), and a(T0)
which is much higher than that of ceramic plate, hence the internal
are the thermal conductivity, electric resistivity, and Seebeck coef-
temperature difference of heat sink is expected to be very small. The
ficient at reference temperature, respectively; A1, A2, B1, B2, C1,
simulations verify the guess by the fact that the internal tempera-
and C2 are the fitted parameters determined by experimental data
ture difference is less than 0.2 K, therefore, we have T top
sink ¼ T max .
as listed in Table 1.
The equivalent heat transfer coefficient can be calculated from
The connector and heat sink are made of copper and their prop-
Eqs. (19) and (20), or:
erties are assumed to be constant due to weak thermoelectric ef-
fect. AlN is selected as material of ceramic plates and its
properties are listed in Table 2 [30]. QH q Hf D2g Dp
heff ¼ ¼ a tf cp;a ð21Þ
AðT max  T 1 Þ 1 þ 12la L2sink
Dg
3. Simplification of heat sink
In the present work, with Dp = 1.3 Pa, heff = 7712 W m2 K1.
For a TEC, it is desirable to achieve a high cooling capacity and a
large coefficient of performance (COP) through selecting a proper 380
applied current under the same heat sink cooling condition. Thus,
for a fair comparison, the heat sink is maintained at the same oper- 370
ation conditions, which means that the equivalent convective heat
360
transfer coefficients heff in Eq. (10) are the same for all transient
simulations. 350
The most common configuration of heat sink used in current
340
applications is of the fin structure with the geometry shown in
T (K)

Fig. 2. The heat sink comprises a series of parallel fins with height 330
Hf, thickness tf, and length Lsink = L2. Each fin is spaced by a gap Dg
and mounted on the heat sink base with an area of Lsink  Wsink. Air 320
is employed as the coolant. The geometric parameters of the heat Δt = 0.1 s
sink and the air properties are shown in Table 3 [31,32]. In order
310 Δt = 0.5 s
Δt = 1 s
to simplify the calculation, the heat sink with real structure is sim- 300
Δt = 2 s
plified to a flat plate with the same area, which was also adopted
290
by Chen et al. [33]. In order to ensure the same heat dissipation 0 200 400 600 800 1000 1200 1400 1600
for the equivalent and real heat sinks, the equivalent convective t (s)
heat transfer coefficient heff between the equivalent heat sink
and the ambient must meet: Fig. 3. Independence examination of time step size.
344 J.-H. Meng et al. / Applied Energy 108 (2013) 340–348

306 load test is to avoid the errors in measuring the magnitude of the
304 (a) I = 0.3 A cooling load and the heat losses with it which could lead to uncer-
tainties in the testing conditions [26]. Fig. 4 compares the measure-
302
ment data regarding the dynamic temperature variations at the cold
300 and hot ends with numerical predictions by the present model at
298 I = 0.3, 0.5 and 1 A, respectively. It is found that the predictions agree
296 well with the experiments, the maximum temperature difference
T (K)

294
between numerical predictions and the measurement data is less
than 1.5 K at I = 0.3 A, 2.0 K at I = 0.5 A, and 4.5 K at I = 1 A, respec-
292
tively, hence, the present transient model can accurately analyze
290 the dynamic characteristics of TECs.
288
Experimental data [35]
286
Numerical predictions
284
5. Simulation cases

0 100 200 300 400 500

Thermal conductivity, electrical conductivity, and Seebeck coef-
t (s)
ficient of semiconductor materials are all temperature-dependent.
315 Our previous study [29] has shown that these properties signifi-
(b) I = 0.5 A TL cantly affect the cooling capacity and the COP of TECs when a large
312
temperature gradient occurs within semiconductor materials. For
309 the n-type and p-type semiconductors of Bi2(Te0.94Se0.06)3 and
306 (Bi0.25Sb0.75)Te3, with a small applied current, there is no marked
difference in the cooling capacity and COP predicted by constant
303
property model and variable property model, however, as the cur-
300
T (K)

rent is larger than 1.2 A, ignoring the property change overesti-

TH
297 mates the cooling capacity and COP. Especially, note that the
294 constant property model significantly overestimates the maximum
working current of TEC, which is defined as the corresponding
291
maximum applied current when the cooling capacity is greater
288 Experimental data [35] than zero.
285 Numerical predictions In the present work, the TEC geometry and the convective heat
282
transfer coefficient of heat sink remain unchanged, with the same
0 100 200 300 400 500 cooling load at the cold end, the response time of temperatures
t (s) (defined as the consuming time for the temperature of the cold
and hot ends reaches the steady-state), the dynamic temperature
380
variations, as well as the steady-state temperatures and tempera-
370 (c) I = 1 A ture difference at the cold and hot ends are closely dependent on
360 the applied current. At different applied currents, the Joule heat
and Thomson heat generated within semiconductors have large
350
difference, which leads to a remarkable difference in temperature
340 distributions of semiconductors. Hence, it can be expected that
T (K)

330 the variable property model could predict the TEC dynamic charac-
320 teristics more accurately when a large temperature gradient occurs
within semiconductors. Therefore, this work focuses on the com-
310
parisons of the dynamic characteristics of the TEC predicted by
300 constant and variable property models under different applied cur-
290 Experimental data [35] rents and different cooling loads at the cold end. The parameters
Numerical predictions used in the simulations are listed in Table 4, where I = 0.6 A,
280
1.5 A, 2.0 A and 2.8 A denote typical small, medium and large ap-
270
0 100 200 300 400 500 plied currents, respectively, and qL is varied from 14,285.7 W m2
t (s) to 171,428.6 W m2 denoting typical low, medium and high
cooling loads at the cold end.
Fig. 4. Comparison of numerical predictions and experimental data: (a) I = 0.3 A;
(b) I = 0.5 ; (c) I = 1.0 A.
5.1. At small applied current

4. Model validation Fig. 5 shows the dynamic characteristics predicted by the con-
stant and variable property models under different cooling loads
The model used uniformly distributed elements in the x, y and z at small applied current of I = 0.6 A. It can be observed that the re-
directions. The grid independence was examined in preliminary sponse time, the temperature difference between the cold and hot
test runs. The independence of the time step size was also exam- ends, and the temperatures at the cold and hot ends predicted by
ined with Dt = 0.1, 0.5, 1, and 2 s, as shown in Fig. 3, and then both models are almost identical. These can be explained by the
Dt = 1 s was selected. fact that at small current the Joule heat and the Thomson heat gen-
Cheng et al.’s experimental data [34] are used to validate the erated inside the semiconductors are very low, leading to a small
present transient model. A zero-cooling-load test is adopted in the temperature gradient within the semiconductors, and hence a
present comparison. The reason for performing the zero-cooling- invisible property change.
 J.-H. Meng et al. / Applied Energy 108 (2013) 340–348 345

Table 4
Cases simulated in the present work.

Cases I (A) qL1 (W m2) qL2 (W m2) qL3 (W m2) qL4 (W m2) qL5 (W m2)
Low current 0.6 0 14285.7 28571.4 42857.1 71428.6
Medium current 1.5 0 57142.9 85714.3 114285.7 142857.1
2.0 0 28571.4 71428.6 142857.1 171428.6
Large current 2.8 0 28571.4 71428.6 142857.1 171428.6

Fig. 5 shows that the cold end temperature decreases but the then increases until reaching the steady-state. This phenomenon
hot end temperature increases with the time increased. However, can be explained as follows: the Peltier cooling occurs at the cold
both gradually approach to the steady-state. With the increase in end of the semiconductors, while the Joule heating and Thomson
the cooling load, the temperatures at the cold and hot ends are effect occurs uniformly across the semiconductors columns. Dur-
both increased, but the temperature rising rate of the cold end is ing the dynamic process, the cooling at the cold end occurs before
much higher than that of the hot end, resulting in a significant the Joule heat and Thomson heat reach this end and a temporary
reduction in the temperature difference between the cold and temperature overshoot at the cold end can be observed. However,
hot ends. This is responsible for that with the same heff and I, only for small qL, there is a large temperature differences between the
a low temperature difference can reduce Fourier’s thermal conduc- cold and hot ends, leading to a enhanced Fourier’s back heat con-
tion, hence the TEC maintains a high cooling capacity. It is noted duction from the hot end to the cold end, hence the undershoot
that, since the temperature difference decreases and the does not occur.
steady-state temperature at the cold end is closer to the initial Fig. 8 shows the electric potential difference and COP (= QL/
temperature of 300 K under high cooling load, the TEC response P = QL/IV, where QL = qLA is the cooling capacity of the TEC, P is
time is significantly reduced as the cooling load increases. the electric power, V is the electric potential difference) predicted
by the constant and variable property models at I = 1.5 A. It can be
seen that due to the negative temperature effect of electric conduc-
5.2. At medium applied current tivity for Bi2(Te0.94Se0.06)3 and (Bi0.25Sb0.75)Te3 materials, the vari-
able property model predicts higher electric potential difference
Fig. 6 shows the dynamic temperature variations at the cold and through the TEC. Therefore, more electric power is needed for the
hot ends predicted by the constant and variable property models at variable property model, leading to a reduced COP.
I = 1.5 A with various cooling loads. With zero cooling load, the dy- Fig. 9 shows the dynamic temperature variations at the cold and
namic temperature variations at the cold end are almost the same hot ends at I = 2.0 A with qL = 0 W m2 and 142,857.1 W m2. It is
with that at the hot end, only a slightly difference occurs when noted that at the larger applied current of I = 2.0 A, there is a
t > 200 s. The two models both predict the temperature overshoot remarkable difference between the constant and variable models
phenomenon at the hot end, for example, the variable property even at zero cooling load due to significantly enhanced Joule heat.
model indicates that the hot end temperature increases rapidly With high cooling load qL = 142857.1 W m2, the temperature at
from the initial stage of 300 K, then decreases slowly after reaching the cold end predicted by the variable property model does not de-
the maximum value of 348.1 K, finally approaches to the steady va- creases but increases continuously to steady-state value of 392.4 K,
lue of 342 K at the time of 1200 s. which is even higher than the temperature at the hot end predicted
The constant and variable property models predict that, at by constant property model. The results again confirm that tem-
I = 1.5 A the temperatures at the cold and hot ends both increase perature-dependent properties have significant effect on dynamic
with the increase in the cooling load, which has been observed at characteristics of the TEC.
small current of 0.6 A (Fig. 5). However, at medium current of
1.5 A, the difference in dynamic temperature variation predicted
by both models occurs and the difference become more remark-
able at high cooling load. The constant property model underesti-
mates the temperatures at the cold and hot ends. In details, the
steady-state temperatures at the hot end and cold end predicted
by the constant property model under qL = 85,714.3 W m2 is
330 Constant property: red line
347.8 K, and 281.6 K, which underestimates by 6.8 K and 6.7 K Variable property: black line TH
compared to the variable property model, respectively. Under
315
qL = 142,857.1 W m2, however, the underestimations increase to
11.0 K and 14.7 K at the hot end and cold end, respectively. Fig. 7 300
shows the temperature distributions within p- and n-type semi-
conductors predicted by both models under qL = 142,857.1 W m2, 285
T (K)

and t = 1000 s. It can be seen that the temperature predicted by the TL

constant property model is lower than that by the variable prop- 270
erty model. According to Eqs. (2) and (3), the reason can be ex-
255
plained. Because the electrical conductivities of Bi2(Te0.94Se0.06)3
: qL= 0 W m-2
and (Bi0.25Sb0.75)Te3 have negative temperature effect, the constant 240
property model underestimates the Joule heat generated within Ο : qL= 28571.4 W m-2
the semiconductors, leading to a lower temperature distribution 225 Δ : qL=71428.6 W m-2
(Fig. 7). 0 200 400 600 800 1000 1200 1400 1600 1800 2000
Note that with high cooling loads of qL = 114,285.7 W m2 and
t (s)
142,857.1 W m2, a undershoot phenomenon of cold end temper-
ature is observed (Fig. 7c and d). The cold end temperature de- Fig. 5. Dynamic temperature variations at the cold and hot ends predicted by
creases from 300 K to the minimum value at the initial stage, constant and variable property models at I = 0.6 A.
346 J.-H. Meng et al. / Applied Energy 108 (2013) 340–348

(a) 350 (b) 360

340 350
330 TH TH
340
320
310 330
300
320 Constant Property
T (K)

T (K)
290 Constant Property
Variable Property Variable Property
280 310
270
300
260 TL
250
TL 290
240
280
230 I = 1.5A qL=0 W m -2 I = 1.5A qL=85714.3 W m-2
220 270
0 200 400 600 800 1000 1200 1400 1600 1800 2000 0 200 400 600 800 1000 1200 1400 1600 1800 2000
t (s) t (s)

(c) 360
(d) 370

TH 360
TH
350
350
340 114285.7
340 Constant Property
330 Variable Property 1428
T (K)

T (K)
Constant Property
330
Variable Property
320
320 TL
310
310
300 TL
300
290 I = 1.5A qL=114285.7 W m-2 I = 1.5A qL=142857.1 W m-2
290
0 200 400 600 800 1000 1200 1400 1600 1800 2000 0 200 400 600 800 1000 1200 1400 1600 1800 2000
t (s) t (s)

Fig. 6. Dynamic temperature variations at the cold and hot ends predicted by constant and variable property models at I = 1.5 A: (a) qL = 0 W m2, (b) qL = 85,714.3 W m2, (c)
qL = 114,285.7 W m2, (d) qL = 142,857.1 W m2.

qL = 0 W m2, while TL = 348.3 K and TH = 438.8 K at

390 I =1.5 A,t=1000 s, qL =142857.1 W m-2
qL = 142,857.1 W m2 when the TEC element approaches the stea-
380 dy-state. However, the variable property model indicates that the
temperatures at the cold and hot ends increase rapidly, and it
370
can never reach the steady-state
360 Using our previous steady-state model [29], the I–QL perfor-
mance curves with fixed temperatures at the cold and hot ends,
350
T (K)

TL = 295.0 K and TH = 430.2 K, as well as TL = 348.3 K and

340 TH = 438.8 K, can be predicted as shown in Fig. 11. The black lines
are predicted by constant property model, and red dash dot lines
330 p-type, constant properties
are predicted by variable property model. Note that the hollow cir-
p-type, variable properties
320 n-type, cosntant properties cle in I–QL curve (constant property model, TL = 295.0 K and
n-type, variable properties
310

0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 0.39
0.18 I=1.5 A,qL=142857.1 W m -2
z (mm)
0.16 0.36
Fig. 7. Temperature distributions within semiconductors predicted by constant and V, constant property
0.14
variable property models at I = 1.5 A. V, variable property 0.33
0.12
0.30
0.10
COP
V (V)

5.3. At large applied current 0.27

0.08
0.24
Fig. 10 shows the dynamic temperature variations at the cold 0.06
and hot ends at I = 2.8 A with qL = 0 W m2 and 142,857.1 W m2, 0.21
0.04 COP, constant property
respectively. At large applied current, the temperature gradient COP, variable property
within the semiconductors is significantly elevated due to the 0.02 0.18
strong Joule heating, resulting in strong effect of material proper- 0.00 0.15
ties. Consequently, dynamic temperature variations at the cold 0 200 400 600 800 1000 1200 1400 1600 1800 2000
and hot ends predicted by constant property model remarkably t (s)
differ from those predicted by variable property model. The con-
Fig. 8. The electric potential difference and COP predicted by constant and variable
stant property model predicts TL = 295.0 K and TH = 430.2 K at
property models at I = 1.5 A.
 J.-H. Meng et al. / Applied Energy 108 (2013) 340–348 347

(a) 405 (b) 450

435
390
TH 420 TH
375
360 405 TL
345 390

330 375 TH
T (K)

T (K)
Constant property
315 Variable property 360
300 345 Constant property
Variable property
285 330
270 315
TL TL
255 300
I = 2.0 A qL= 0 W m 2 I = 2A qL= 28571.4 W m 2
240 285
0 200 400 600 800 1000 1200 1400 1600 1800 2000 0 200 400 600 800 1000 1200 1400 1600 1800 2000
t (s) t (s)

Fig. 9. Dynamic temperature variations at the cold and hot ends predicted by constant and variable property models at I = 2.0 A: (a) qL = 0 W m2, (b) qL = 142,857.1 W m2.

(a) 500 (b) 500 2

480 I = 2.8A q L= 0 W m
2 Constant property 480 I =2.8A, qL=142857.1 W m Constant property
Varivable property Varivable property
460 TH TL 460 TH TL
440 440
420 TH TH
420
400
400
T (K)

T (K)

380
380
360
360
340
320 340 TL
300 320
280 TL 300
260 280
0 200 400 600 800 1000 1200 1400 1600 1800 2000 0 200 400 600 800 1000 1200 1400 1600 1800 2000
t (s) t (s)

Fig. 10. Dynamic temperature variations at the cold and hot ends predicted by constant and variable property models at I = 2.8 A: (a) qL = 0 W m2, (b) qL = 142,857.1 W m2.

constant property
working current ranges from 0.49 A to 2.40 for fixed TL = 295.0 K
0.10 TL=295.0 K and TH = 430.2 K, and from 1.09 A to 1.94 A for fixed TL = 295 K
TH=430.2 K and TH = 430.2 K. Apparently, for the transient cases, applied current
I = 2.8 A is higher than the maximum working current. Under this
0.08 variable property
TL=295.0 K circumstance, the heat transferred from the hot end to the cold
TH=430.2 K end by Fourier’s heat conduction is still larger than the heat ad-
0.06 sorbed at the cold end by Peltier effect even at zero cooling load,
QL (W)

thus, the temperatures at the cold and hot ends both increase result-
ing the failure of TEC refrigeration.
0.04
variable property
TL=348.3 K constant property
TH=438.8 K TL=348.3 K
0.02 6. Conclusions
TH=438.8 K

This work develops a complete three-dimensional transient TEC

0.00
0.0 0.4 0.8 1.2 1.6 2.0 2.4 2.8 3.2 3.6 4.0 model, which introduces the coupling of the heat transfer and elec-
I (A) tric conduction within semiconductors, and the transient TEC mod-
el developed in this paper has not been found in the existing
Fig. 11. Cooling capacity of TEC at different temperature difference between the literatures. The developed model is used to investigate the dy-
cold and hot ends predicted by constant and variable property models. namic characteristics of TECs. A miniature TEC element with a
finned structure heat sink is investigated and Bi2(Te0.94Se0.06)3
and (Bi0.25Sb0.75)Te3 are selected as the n-type and p-type thermo-
TH = 430.2 K) denotes the steady-state point with I = 2.8 A and electric materials, respectively. The present model predicts that
QL = 0.1 W (qL = QL/A = 142,857.1 W m2), the solid circle in I–QL dynamic temperature variations at cold and hot ends under zero
curve (constant property model, TL = 348.3 K and TH = 438.8 K) cooling load agree well with experimental results.
denotes the steady-state point with I = 2.8 A, QL = 0 W. The result Three types of operating conditions including small, medium,
indicates that the transient model developed in this work can and large applied currents are adopted to examine the effect of
predict well the dynamic response characteristics and the final stea- material properties on the dynamic characteristics of the TEC.
dy-state. The variable property model predicts that the effective The results show that at small currents the constant and variable
348 J.-H. Meng et al. / Applied Energy 108 (2013) 340–348

property models predict almost the same dynamic characteristics; [11] Lee KH, Kim OJ. Analysis on the cooling performance of the thermoelectric
micro-cooler. Int J Heat Mass Transfer 2007;50:1982–92.
at medium currents the predictions of two models differ, and the
[12] Yu J, Wang B. Enhancing the maximum coefficient of performance of
difference is more apparent at high cooling load; at large currents thermoelectric cooling modules suing internally cascaded thermoelectric
the variable property model predict that the temperatures at cold couples. Int J Refriger 2009;32:32–9.
and hot ends increase continuously and never reach the steady- [13] Chen WH, Liao CY, Huang CI. A numerical study on the performance of
miniature thermoelectric cooler affected by Thomson effect. Appl Energy
state, oppositely, the constant property model predict that the 2012;89:464–73.
TEC still supply refrigeration. The present simulations indicate that [14] Small-scale energy harvesting. In: Lallart M. Three dimensional TCAD
the variable property model is more appropriate to accurate pre- simulation of a thermoelectric module suitable for use in a thermoelectric
energy harvesting system, Could C, Shammas N, InTech; 2012. [Chapter 2]
dict TEC transient performance, especially when the TEC operates [15] Gou XL, Xiao H, Yang SW. Modeling, experimental study and optimization on
at medium and high applied currents. low-temperature waste heat thermoelectric generator system. Appl Energy
At a specific applied current, the temperatures at the cold and 2010;87:3131–6.
[16] Kim SH. Analysis and modeling of effective temperature differences and
hot ends both increase with the elevated cooling load, but the tem- electrical parameters of thermoelectric generators. Appl Energy
perature rising rate of the cold end is much higher than that of the 2013;102:1458–63.
hot end, resulting in a significant reduction in the temperature dif- [17] Suter C, Jovanovic ZR, Steinfeld A. A 1 kWe thermoelectric stack for geothermal
power generation – modeling and geometrical optimization. Appl Energy
ference between the cold and hot ends. This phenomenon can be 2012;99:379–85.
attributed to that small temperature difference reduces Fourier’s [18] van Sark WGJHM. Feasibility of photovoltaic – thermoelectric hybrid modules.
back heat transfer and hence maintains the high cooling capacity Appl Energy 2011;88:2785–90.
[19] Qiu K, Hayden ACS. Integrated thermoelectric and organic Rankine cycles for
of TECs. The response time of the TEC varies from several hundreds
micro-CHP systems. Appl Energy 2012;97:667–72.
of seconds to several thousands of seconds, which is dependent on [20] Patyk A. Thermoelectric generators for efficiency improvement of power
its operating conditions. generation by motor generators – environmental and economic perspectives.
Appl Energy 2013;102:1448–57.
[21] He W, Su YH, Riffat SB, Hou JX, Ji J. Parametrical analysis of the design and
Acknowledgement performance of a solar heat pipe thermoelectric generator unit. Appl Energy
2011;88:5083–9.
[22] Hsu CT, Huang GY, Chu HS, Yu B, Yao DJ. Experiments and simulations on low-
This study was supported by the National Natural Science Foun- temperature waste heat harvesting system by thermoelectric power
dation of China (No. 51276060), by National Basic Research Pro- generators. Appl Energy 2011;88:1291–7.
[23] Alata M, Al-Nimr MA, Naji M. Transient behavior of a thermoelectric device
gram of China (No. 2009CB219803), by Program for New Century
under the hyperbolic heat conduction model. Int J Thermophys
Excellent Talents in University (No. NCET-11-0635), and by the 2003;26:1753–68.
Fundamental Research Funds for the Central Universities (No. [24] Montecucco A, Buckle JR, Knox AR. Solution to the 1-D unsteady heat
11ZG01). conduction equation with internal Joule heat generation for thermoelectric
devices. Appl Therm Eng 2012;35:177–84.
[25] Taylor PJ, Jesser WA, Rosi FD, Derzko Z. A model for the non-steady-state
temperature behaviour of thermoelectric cooling semiconductor devices.
References Semiconductor Sci Technol 1997;12:443–7.
[26] Cheng CH, Huang SY, Cheng TC. A three-dimensional theoretical model for
[1] Riffat SB, Ma X. Thermoelectrics: a review of present and potential predicting transient thermal behavior of thermoelectric coolers. Int J Heat
applications. Appl Therm Eng 2003;23:913–35. Mass Transfer 2010;53:2001–11.
[2] Rowe DM, editor. Thermoelectrics handbook: macro to nano. Boca Raton: CRC [27] Thonhauser T, Mahan GD, Zikatanov L, Roe J. Improved supercooling in
Press; 2006. transient thermoelectrics. Appl Phys Lett 2004;85:3247–9.
[3] Minnich AJ, Dresselhaus MS, Ren ZF, Chen G. Bulk nanostructured [28] Yang RG, Chen G, Ravi Kumar A, JeffreySnyder G, Jean-Pierre Fleurial. Transient
thermoelectric materials: current research and future prospects. Energy cooling of thermoelectric coolers and its applications for microdevices. Energy
Environ Sci 2009;2:466–79. Conv Manage 2005;46:1407–21.
[4] Huang BJ, Duang CL. System dynamic model and temperature control of a [29] Wang XD, Huang YX, Cheng CH, Lin DTW. A three-dimensional numerical
thermoelectric cooler. Int J Refriger 2000;23:197–207. modeling of thermoelectric device with consideration of coupling of
[5] Völklein F, Min G, Rove DM. Modeling of microelectromechanical temperature field and electric potential field. Energy 2012;47:488–97.
thermoelectric cooler. Sensors Actuators 1999;75:95–101. [30] Qin ML, Qu XH, Lin JL, Xiao PA, Zhu BJ, Tang CF. Progress in research and
[6] Huang BJ, Chin CJ, Duang CL. A design method of thermoelectric cooler. Int J development of aluminum nitride (AlN) ceramics. Rare Metal Mater Eng
Refriger 2000;23:208–18. 2002;31:8–12.
[7] Xuan XC, Ng KC, Yap C, Chua HT. A general model for studying effects of [31] Mereu S, Sciubba E, Bejan A. The optimal cooling of a stack of heat generating
interface layers on thermoelectric devices performance. Int J Heat Mass boards with fixed pressure drop, flowrate or pumping power. Int J Heat Mass
Transfer 2002;45:5159–70. Transfer 1993;36:3677–86.
[8] Huang MJ, Yen RH, Wang AB. The influence of the Thomson effect on the [32] Incropera FR, Dewitt DP, Bergman TL, Lavine AS. Fundamentals of heat and
performance of a thermoelectric cooler. Int J Heat Mass Transfer mass transfer. 6th ed. Hoboken: John Wiley & Sons; 2007.
2005;48:413–8. [33] Wang CC, Hung CI, Chen WH. Design of heat sink for improving the
[9] Cheng YH, Shih C. Maximizing the cooling capacity and COP of two-stage performance of thermoelectric generator using two-stage optimization.
thermoelectric coolers through genetic algorithm. Appl Therm Eng Energy 2012;39:236–45.
2006;26:937–47. [34] Chen CH, Huang SY. Development of a non-uniform-current model for
[10] Pan Y, Lin B, Chen J. Performance analysis and parametric optimal design of an predicting transient thermal behavior of thermoelectric coolers. Appl Energy
irreversible multi-couple thermoelectric refrigerator under various operating 2012;100:326–35.
conditions. Appl Energy 2007;84:882–92.