Calculation Methods for Thermoelectric Generator Performance

Abstract

This chapter aims to build one-dimensional thermoelectric model for device-level thermoelectric generator (TEG) performance calculation and prediction under steady heat transfer. Model concept takes into account Seebeck, Peltier, Thomson effects, and Joule conduction heat. Thermal resistances between heat source, heat sink, and thermocouple are also considered. Then, model is simplified to analyze influences of basic thermal and electrical parameters on TEG performance, when Thomson effect is neglected. At last, an experimental setup is introduced to gauge the output power and validate the model. Meantime, TEG simulation by software ANSYS is introduced briefly.

1. Introduction

Output power P_out and energy conversion efficiency η are the primary parameters to characterize TEG performance. They are intensively influenced by such factors as temperature of heat source and sink, thermoelectric materials physical properties, thermocouple geometries, thermal and electrical contact properties, and load factor. Therefore, it is necessary to build physical model formulating these factors concisely, to conduct realistic TEG design. At present, many significant works have been undertaken for modeling device-level TEG precisely [1–3]. In addition, comprehensive three-dimensional (3D) thermoelectric model has been successfully developed in software ANSYS [4]. In Refs. [5–7], quasi-one-dimensional thermoelectric model is established, where Thomson effect and thermal resistances between thermocouple and heat source, heat sink are neglected. In Ref. [8], improved one-dimensional model including Thomson coefficient and thermal resistances is used to analyze the matched load, the limit of energy conversion efficiency, and the influence of Peltier effect. It shows, that expression of matched load contains not only the inner electrical resistance of TEG, but also the terms resulting from Peltier and Joule effects. In Ref. [9], one-dimensional model to analyze the influence of Thomson heat is built and experimentally validated.

In this chapter, Seebeck, Peltier, Thomson effect, and Joule conduction heat are formulated in thermoelectric generation module model. By model simplification, analytical expressions of output power and energy efficiency are introduced. Essential factors for enhancing the output power are extracted. Then, an experimental setup is built to measure the output power and validate the model. And TEG simulation by software ANSYS is presented.

2. Thermoelectric model for device-level TEG

2.1. TEG cell structure

TEG cell consisting of thermocouple is shown in Figure 1, where basic thermoelectric effects including Peltier and Joule heat and a circuit with load R_L are included. The p and n thermoelements are cuboids of the same thickness and bridged by an electrode in series. Practical devices usually make use of thermoelectric modules containing a number of TEG cells connected electrically in series and thermally in parallel. Cross-sectional area and thickness of thermocouple are marked as A and l. Subscripts ‘n’ and ‘p’ are used to discriminate conductivity type of thermoelements. Temperature of heat source and heat sink is T₁ and T₀, and that of hot and cold side of thermocouple is T_h and T_c. ∆T_g = T_h − T_c is temperature difference on thermoelements, and ∆T = T₁ − T₀ is the one of heat source and heat sink.

Figure 1.

Structure and circuit sketch of TEG cell.

There are Joule heat flowing out and Peltier heat flowing in at hot end of thermoelements, and at cold end, Peltier heat flows out and Joule heat flows out. In addition, there is thermal resistance R_th,h and R_th,c between thermoelements, heat source and heat sink. Heat flow q_h passes from heat source to hot side of thermocouple and the counterpart q_c outflows from cold side of thermocouple to heat sink.

2.2. Basic model

It is assumed, that thermoelements are physically homogeneous and insulated from the surroundings both electrically and thermally, except at junction-reservoir contacts [8–9]. Variable x is defined as location in the thickness direction of thermoelements. According to nonequilibrium thermodynamics under steady heat transfer, energy conservative equations of temperature distributions T_n(x) and T_p(x) are:

$\{\begin{array}{l}{K}_{\text{n}}{l}_{\text{n}}\frac{d^2{T}_{\text{n}}(x)}{d{x}^2}-{\tau }_{\text{n}}I\frac{d{T}_{\text{n}}(x)}{dx}+\frac{R_{\text{n}}{I}^2}{l_{\text{n}}}=0\\ K_{\text{p}}{l}_{\text{p}}\frac{d^2{T}_{\text{p}}(x)}{d{x}^2}+{\tau }_{\text{p}}I\frac{d{T}_{\text{p}}(x)}{dx}+\frac{R_{\text{p}}{I}^2}{l_{\text{p}}}=0\end{array}.$E1

Three terms in the above equations represent thermal conduction, Thomson and Joule heat. K, R, and τ are thermal conductance, electrical resistance and Thomson coefficient (V∙K⁻¹), respectively. Relationship of K, R and A, l is $K=\frac{\lambda A}{l}$ and $R=\frac{\rho l}{A}$, where λ and ρ are thermal conductivity and electrical resistivity of thermoelectric materials. To solve Eq. (1) analytically, material parameters K, R, and τ are considered to be constant. The boundary conditions of Eq. (1) are:

$T_{\text{n}}\left(0\right)=T_{\text{p}}\left(0\right)=T_{\text{c}},$E2

$T_{\text{n}}\left(l_{\text{n}}\right)=T_{\text{p}}\left(l_{\text{p}}\right)=T_{\text{h}}.$E3

Electrical current I is determined by formula:

$I=\frac{U_0}{R_{\text{L}}+R_{\text{g}}},$E4

$U_0={\displaystyle {\int }_{T_{\text{c}}}^{T_{\text{h}}}\alpha (T)\text{d}T},$E5

where U₀ is the voltage of thermocouple, R_g is the electrical resistance of TEG cell, which contains resistance of thermocouple and contact resistance, and α(T) = α_p(T) − α_n(T) is Seebeck coefficient (V∙K⁻¹) of thermocouple.

In practice, temperature of heat source T₁ and heat sink T₀ can be measured and determined. To acquire T_h and T_c, relationship of T₁, T₀ and T_h, T_c is necessary. That is:

$\{\begin{array}{l}{T}_{\text{1}}-T_{\text{h}}=R_{\text{t}}{}_{\text{h,h}}{q}_{\text{h}}\\ T_{\text{c}}-T_{\text{0}}=R_{\text{t}}{}_{\text{h,c}}{q}_{\text{c}}\end{array}.$E6

In Eq. (6), heat flows q_h and q_c are:

$q_{\text{h}}=K_{\text{n}}{l}_{\text{n}}\frac{\text{d}{T}_{\text{n}}\left(x\right)}{\text{d}x}|{}_{x=l_{\text{n}}}+K_{\text{p}}{l}_{\text{p}}\frac{\text{d}{T}_{\text{p}}\left(x\right)}{\text{d}x}|{}_{x=l_p}+\alpha (T_{\text{h}})T_{\text{h}}I-I^2{R}_{\text{ch}},$E7

$q_{\text{c}}=K_{\text{n}}{l}_{\text{n}}\frac{\text{d}{T}_{\text{n}}\left(x\right)}{\text{d}x}|{}_{x=0}+K_{\text{p}}{l}_{\text{p}}\frac{\text{d}{T}_{\text{p}}\left(x\right)}{\text{d}x}|{}_{x=0}+\alpha (T_{\text{c}})T_{\text{c}}I+I^2{R}_{\text{cc}},$E8

wherein R_ch and R_cc are contact electrical resistances at hot and cold side of the thermocouple. Thermal conduction heat, Peltier heat (the third term), and contact Joule heat are within Eqs. (7) and (8). By solving Eq. (1) with Eqs. (2)–(5), T_n(x) and T_p(x) only relating to T_h, T_c, and R_L can be obtained. And flows q_h and q_c can be formulated with T_h, T_c, and R_L in Eqs. (7) and (8). Then, T_h and T_c can be determined for a given R_L by solving Eq. (6) numerically, which is presented in detail [9].

Finally, the output power P_out and energy conversion efficiency η are calculated by the basic equations of thermoelectricity:

$P_{\text{out}}=q_{\text{h}}-q_{\text{c}}=I^2{R}_{\text{L}},$E9

$\eta =\frac{P}{q_{\text{h}}}.$E10

When neglecting Thomson heat, the problem will be much simplified. By solving Eqs. (1)–(8) with τ = 0, an cubic equation about ΔT_g can be yielded as:

$a_1\Delta T_g^3+b_1\Delta T_g^2+c_1\Delta T_g+1=0.$E11

where:

$a_1=\frac{{\overline{\alpha }}^4{R}_{\text{L}}{R}_{\text{th,c}}{R}_{\text{th,h}}}{{\left(R_{\text{g}}+R_{\text{L}}\right)}^3\Delta T}$E26

$b_1=-\frac{{\overline{\alpha }}^2{R}_{\text{th,}}{}_{\text{c}}{R}_{\text{th,}}{}_{\text{h}}}{{\left(R_{\text{g}}+R_{\text{L}}\right)}^2\Delta T}\left[R_{\text{g}}\left(\frac{\epsilon }{R_{\text{th,}}{}_{\text{c}}}+\frac{\epsilon -1}{R_{\text{th,}}{}_{\text{h}}}\right)+R_{\text{L}}\left(\frac{1}{R_{\text{th,}}{}_{\text{c}}}-\frac{1}{R_{\text{th,}}{}_{\text{h}}}\right)\right],$E27

$c_1=-\frac{R_{\text{th,c}}{R}_{\text{th,h}}}{\Delta T}\left[\left(\frac{T_1}{R_{\text{th,h}}}+\frac{T_0}{R_{\text{th,c}}}\right)\frac{{\alpha }^2}{R_{\text{g}}+R_{\text{L}}}+\frac{K_{\text{g}}}{R_{\text{th,h}}}+\frac{K_{\text{g}}}{R_{\text{th,c}}}+\frac{1}{R_{\text{th,c}}{R}_{\text{th,h}}}\right],$E28

$\epsilon =\frac{0.5R_{\text{0}}+R_{\text{cc}}}{R_{\text{g}}},$E29

and Seebeck coefficient α becomes a constant. Equation (11) suits thermoelectric module consisting of m thermocouples, as well, where α and K_g are m times of those of a single thermocouple, but R_th,h and R_th,c are exactly on the contrary.

Generally, c₁ is far larger than a₁ and b₁ in absolute value. Mainly, because of practical module, Seebeck coefficient α has a very small value of about 10⁻² V∙K⁻¹, which is much less than unity. For example, taking module TEG-127-150-9 in Ref. [8], α = 0.05 V∙K⁻¹, R_g = 3.4 Ohm, R_L = 4 Ohm, R_th,c = 6 K∙W⁻¹, R_th,h = 0.1 K∙W⁻¹, K_g = 2.907 W∙K⁻¹, ε ≈ 0.5, T₀ = 297 K, and T₁ = 323 K, calculation result is a₁ ≈ 1.423 × 10⁻⁹ K⁻³, b₁ ≈ 5.905 × 10⁻⁵ K⁻², and c₁ ≈ −0.7461 K⁻¹. So, the terms with ΔT_g order higher than unity can be neglected. At last, here is:

$\Delta T_{\text{g}}=\frac{\Delta T}{1+R_{\text{th,c}}{K}_{\text{g}}+R_{\text{th,h}}{K}_{\text{g}}+{\alpha }^2\frac{\left(R_{\text{th,c}}{T}_1+R_{\text{th,h}}{T}_0\right)}{R_{\text{g}}+R_{\text{L}}}}.$E12

It can be seen, that ΔT_g is influenced not only by thermal resistances R_th,c and R_th,h, but also by Peltier effect, which is presented in the last term of the denominator and functions to decrease ΔT_g. Because it is tantamount to accelerate heat conduction in thermocouple, Peltier heat flows in and out on two sides of thermocouple. By combing Eqs. (12) and (4), (5), (9), the output power P_out is:

$P_{\text{out}}=\frac{{\alpha }^2\Delta T_{\text{g}}^2{R}_{\text{L}}}{{\left(R_{\text{g}}+R_{\text{L}}\right)}^2}=\frac{{\alpha }^2\Delta T^2{R}_{\text{L}}}{{\left(1+R_{\text{th,}}{}_{\text{c}}{K}_{\text{g}}+R_{\text{th,}}{}_{\text{h}}{K}_{\text{g}}\right)}^2{\left[R_{\text{g}}+R_{\text{L}}+\frac{{\alpha }^2\left(R_{\text{th,}}{}_{\text{c}}{T}_1+R_{\text{th,}}{}_{\text{h}}{T}_0\right)}{1+R_{\text{th,}}{}_{\text{c}}{K}_{\text{g}}+R_{\text{th,}}{}_{\text{h}}{K}_{\text{g}}}\right]}^2}.$E13

R_L, T₁, T₀, R_th,c, R_th,h, α, R_g, and K_g directly affect P_out. In those parameters, α, R_g and K_g are TEG internal factors, and R_L, T₁, and T₀ are the external ones, and R_th,c and R_th,h originate from both the internal and external. From the form of Eq. (13), it is obvious, that reducing T₁, T₀, R_th,c, and R_th,h can increase P_out, if ΔT is constant, owing to influence of Peltier effect on ΔT_g. On the other hand, P_out has a maximum along with R_g and K_g.

2.3. Matched load, output power and energy efficiency

First of all, influence of R_L on P_out is analyzed. In Eq. (13), P_out reaches maximum, when R_L is:

$R_L=R_{\text{g}}+\frac{{\alpha }^2\left(R_{\text{th,c}}{T}_1+R_{\text{th,}}{}_{\text{h}}{T}_0\right)}{1+R_{\text{th,c}}{K}_{\text{g}}+R_{\text{th,}}{}_{\text{h}}{K}_{\text{g}}},$E14

which is the matched load and marked as R_L,m. Indeed, R_L,m is slightly larger than R_g due to the very small value of α². It means, that existence of Peltier effect increases irreversible heat in thermoelectric module. And reducing T₁ and T₀ helps to cut down this irreversible heat. When K_g→+∞, R_L,m is equal to R_g, since at this moment heat conduction in thermocouple runs under infinitesimal temperature difference and the irreversibility of heat transfer disappears. However, this irreversibility exists with finite K_g, leading to heat loss in thermocouple, that is equivalent to increase in internal resistance. Define R_L/R_g as the load factor s_L. So, s_L is:

$s_{\text{L}}=1+\frac{Z{K}_{\text{g}}\left(R_{\text{th,c}}{T}_1+R_{\text{th,}}{}_{\text{h}}{T}_0\right)}{1+R_{\text{th,c}}{K}_{\text{g}}+R_{\text{th,}}{}_{\text{h}}{K}_{\text{g}}},$E15

when R_L is equal to matched load and $Z=\frac{{\alpha }^2}{K_{\text{g}}{R}_{\text{g}}}$ is the figure of merit. For thermoelectric module, the output power is:

$P_{\text{out}}=\frac{m^2{\alpha }^2\Delta T^2{R}_{\text{L}}}{{\left(1+R_{\text{th,c}}{K}_{\text{g}}+R_{\text{th,}}{}_{\text{h}}{K}_{\text{g}}\right)}^2{\left\{R_{\text{L}}+m\left[R_{\text{g}}+\frac{{\alpha }^2(R_{\text{th,c}}{T}_1+R_{\text{th,}}{}_{\text{h}}{T}_0)}{1+R_{\text{th,c}}{K}_{\text{g}}+R_{\text{th,}}{}_{\text{h}}{K}_{\text{g}}}\right]\right\}}^2},$E16

where m is the number of thermocouples. And the corresponding matched load is $R_{\text{L,m}}=m\left[R_{\text{g}}+\frac{{\alpha }^2(R_{\text{th,c}}{T}_1+R_{\text{th,}}{}_{\text{h}}{T}_0)}{1+R_{\text{th,c}}{K}_{\text{g}}+R_{\text{th,}}{}_{\text{h}}{K}_{\text{g}}}\right].$

As for energy efficiency η, by Eq. (7), which can be expressed as function of ΔT_g and Eqs. (12) and (13), it is:

$\eta =\frac{P}{q_{\text{h}}}=\frac{{\alpha }^2\Delta T^2{R}_{\text{L}}}{c_2^2{\left(R_{\text{L}}+R_g\right)}^2\left\{\frac{{\alpha }^2\Delta T\left[T_{\text{h}}\left(R_{\text{L}}+R_g\right)-\frac{\epsilon \Delta T{R}_g}{c_2}\right]}{c_2{\left(R_{\text{L}}+R_g\right)}^2}+\frac{\Delta T{K}_{\text{g}}}{c_2}\right\}},$E17

where

$c_2=1+R_{\text{th,c}}{K}_{\text{g}}+R_{\text{th,h}}{K}_{\text{g}}+\frac{{\alpha }^2\left(R_{\text{th,c}}{T}_1+R_{\text{th,h}}{T}_0\right)}{R_{\text{L}}+R_g}.$E30

By solving Eq. (17) about the partial derivative of R_L, it can be obtained, that when load factor s_L is:

$s_{\text{L}}=\sqrt{1+Z{T}_{\text{h}}+\frac{Z\epsilon \Delta T+Z{K}_{\text{g}}\left(1+Z{T}_{\text{h}}\right)\left(R_{\text{th,c}}{T}_1+R_{\text{th,h}}{T}_0\right)}{1+R_{\text{th,c}}{K}_{\text{g}}+R_{\text{th,h}}{K}_{\text{g}}}},$E18

then η reaches maximum. Equation (18) is downright different from Eq. (15) in the expressions, so achieving maximum of output power and energy efficiency simultaneously is impossible. Actually, the corresponding load factor of the former is smaller than that of the latter. When the ideal state is considered (R_th,c = R_th,h = 0), s_L = 1 is for the former and $s_{\text{L}}=\sqrt{1+Z{T}_{\text{h}}+Z\epsilon \Delta T}$, which is larger than 1, is for the latter.

2.4. Influence of K_g on TEG performance

K_g is important internal factor that influences the output performance in TEG. When matched load is reached, the corresponding output power P_out,m is:

$P_{\text{out,m}}=\frac{Z{K}_{\text{g}}\Delta T^2}{4{\left(1+R_{\text{th,}}{}_{\text{c}}{K}_{\text{g}}+R_{\text{th,}}{}_{\text{h}}{K}_{\text{g}}\right)}^2\left[1+\frac{Z{K}_{\text{g}}\left(R_{\text{th,}}{}_c{T}_1+R_{\text{th,}}{}_{\text{h}}{T}_0\right)}{1+R_{\text{th,}}{}_{\text{c}}{K}_{\text{g}}+R_{\text{th,}}{}_{\text{h}}{K}_{\text{g}}}\right]}.$E19

For a common thermoelectric module, thermoelements have the same size, l_n = l_p and A_n = A_p, so the figure of merit Z is not related to their size, but material physical parameters. From Eq. (19), we can see, that increase in Z will enhance the output power. By solving Eq. (19) regarding the partial derivative of K_g, when:

$K_{\text{g}}=\frac{1}{\sqrt{{\left(R_{\text{th,}}{}_{\text{c}}+R_{\text{th,}}{}_{\text{h}}\right)}^2+Z\left(R_{\text{th,}}{}_{\text{c}}{T}_1+R_{\text{th,}}{}_{\text{h}}{T}_0\right)\left(R_{\text{th,}}{}_{\text{c}}+R_{\text{th,}}{}_{\text{h}}\right)}}=\left({\lambda }_{\text{p}}+{\lambda }_{\text{n}}\right)\frac{A{}_{\text{e}}}{l_{\text{e}}}$E20

then P_out,m reaches maximum, where l_e and A_e are thickness and cross-sectional area of thermoelements. Since R_th,c and R_th,h are related to A_e, but not to l_e, there is an optimal l_e to maximize P_out,m:

$P_{\text{out,m}}=\frac{Z\Delta T^2}{4c_3{\left(1+\frac{R_{\text{th,}}{}_{\text{c}}}{c_3}+\frac{R_{\text{th,}}{}_{\text{h}}}{c_3}\right)}^2\left[1+\frac{Z\left(R_{\text{th,}}{}_{\text{c}}{T}_1+R_{\text{th,}}{}_{\text{h}}{T}_0\right)}{c_3\left(1+\frac{R_{\text{th,}}{}_{\text{c}}}{c_3}+\frac{R_{\text{th,}}{}_{\text{h}}}{c_3}\right)}\right]},$E21

and

$c_3=\sqrt{{\left(R_{\text{th,c}}+R_{\text{th,}}{}_{\text{h}}\right)}^2+Z\left(R_{\text{th,c}}^2{T}_1+R_{\text{th,h}}^2{T}_0+R_{\text{th,c}}{R}_{\text{th,}}{}_{\text{h}}{T}_1+R_{\text{th,c}}{R}_{\text{th,}}{}_{\text{h}}{T}_0\right)}.$E31

2.5. Influence of Peltier effect on TEG performance

When Peltier effect is neglected, the relation of ΔT_g and ΔT is:

$\text{Δ}{T}_{\text{g}}=\frac{R_{\text{th,}}{}_g}{R_{\text{th,}}{}_{\text{g}}+R_{\text{th,}}{}_{\text{h}}+R_{\text{th,}}{}_{\text{c}}}\text{Δ}T,$E22

and the corresponding output power with matched load R_L = R_g and constant Seebeck coefficient is:

$P_{\text{out,m}}=\frac{{\left(\alpha \text{Δ}T\right)}^2}{4R_{\text{g}}\left(1+R_{\text{th,h}}{K}_{\text{g}}+R_{\text{th,c}}{K}_{\text{g}}\right)}.$E23

Meantime, the output power Eq. (13) is for the condition without Peltier effect, and the ratio of Eq. (23) and Eq. (19), η^Pelt, reflects influence degree of Peltier effect on the output power:

${\eta }_{\text{Pelt}}={\left(1+\frac{Z{T}_0}{2}\frac{R_{\text{th,}}{}_{\text{c}}+R_{\text{th,}}{}_{\text{h}}}{R_{\text{th,}}{}_{\text{g}}+R_{\text{th,}}{}_{\text{h}}+R_{\text{th,}}{}_{\text{c}}}+\frac{Z}{2}\frac{R_{\text{th,}}{}_{\text{c}}\Delta T}{R_{\text{th,}}{}_{\text{g}}+R_{\text{th,}}{}_{\text{h}}+R_{\text{th,}}{}_{\text{c}}}\right)}^2.$E24

It is known, that when R_th,g ≪ R_th,c + R_th,h, η_Pelt is approximately (1 + 0.5ZT₀ + 0.25ZΔT)² with R_th,c ≈ R_th,h, and even with ΔT → 0, the output power calculated without Peltier effect is more than the output power considering Peltier effect, by over 120% for a common Bi₂Te₃-based module with ZT ≈ 1. That means, the influence of Peltier effect must be considered. Similar status is obtained, where the difference is more than 50%, when R_th,g ≈ R_th,c + R_th,h. On the contrary, when R_th,g ≫ R_th,c + R_th,h, η_Pelt is approximately equal to 1, which means the influence of Peltier effect is negligible. Hence, the smaller the thermal resistance of thermocouple R_th,g, the stronger is the influence of Peltier effect.

Eventually, basic factors for enhancing TEG output power are summarized as:

1. Enhancing ZT of thermoelectric materials and ΔT, decreasing R_th,c and R_th,h.

2. When ΔT is fixed, lower T₀ and T₁ can reduce irreversible heat to elevate output power.

3. Matched load is a little larger than the inner electrical resistance of TEG.

4. There exists an optimal thermocouple thickness to maximize output power.

3. Test validation

3.1. Materials property

P-type and n-type Bi₂Te₃-based materials are, respectively, Bi_0.5Sb_1.5Te₃ and Bi₂Te_2.85Se_0.15, which are prepared by mechanical alloy + spark plasma sintering method. Seebeck coefficient and resistivity of the materials are tested by HGTE-II thermoelectric material performance test system (Chinese patent no. ZL200510018806.4) with test temperature up to 1073 K, relative error of not more than 6%. Thermal conductivity of the materials is measured by laser perturbation method (Type TC-7000 of ULVAC RIKO^®). As shown in Table 1, parameters are obtained by polynomial fitting of the experimental data. In temperature range 273 K < T < 493 K, Seebeck coefficient value α is between 170 × 10⁻⁶ V∙K⁻¹ and 220 × 10⁻⁶ V∙K⁻¹, decreasing with rising temperature. Electrical resistivity ρ is (8.3–20.0) × 10⁻⁶ Ohm∙m and thermal conductivity λ is 1.4–2.1 W∙m⁻¹∙K⁻¹, which both show obvious increase with temperature rise.

Material	Parameters	Values
p-Bi0.5Sb1.5Te3	α/V∙K^−1	−1.791 × 10^−11T^3 + 1.763 × 10^−8T^2 − 5.714 × 10^−6T + 8.304 × 10^−4
	ρ/Ohm∙m	−7.929 × 10^−13T^3 + 7.992 × 10^−10T^2 −1.947 × 10^−7T + 1.728 × 10^−5
	ƛ/W∙m^−1∙K^−1	3.342 × 10^−5T^2 − 2.24 × 10^−2T + 5.118
n-Bi2Te2.85Se0.15	α/V∙K^−1	1.321 × 10^−11T^3 − 1.383 × 10^−8T^2 + 4.81 × 10^−6T − 7.774 × 10^−4
	ρ/Ohm∙m	−7.618 × 10^−13T^3 + 8.098 × 10^−10T^2 − 2.537 × 10^−7T + 3.207 × 10^−5
	ƛ/W∙m^−1∙K^−1	3.264 × 10^−5T^2 − 2.228 × 10^−2T + 5.302

Table 1.

Physical parameters of Bi₂Te₃-based materials (273 K < T < 493 K).

In practice, it is difficult to measure thermal resistances R_th,c and R_th,h, and contact electrical resistances r_cc and r_ch. Their values are determined according to empirical formulas. Contact electrical resistivity ρ_c (Ohm∙m²) at leg-strap junctions and thermal conductivity λ_c (W m⁻¹ K⁻¹) of thermal conductive layer (≈1.2 mm thick) are according to the empirical formulas given by Rowe et al. [5]:

$\frac{2{\rho }_c}{\rho }=0.1\text{mm,}\frac{\lambda }{{\lambda }_c}=\mathrm{0.2.}$E25

Here ρ and λ are electrical resistivity and thermal conductivity of thermoelements, respectively. In our calculation, the mean values of ρ and λ over the temperature range are taken as references. As ρ and λ vary with temperature, values of ρ_c and λ _c are also different as the temperature varies. Experiments under four temperature conditions are carried out and the corresponding parameter values are shown in Table 2.

Temperatures, °C Parameters	T0 = 23; T1 = 81	T0 = 23; T1 = 111	T0 = 27; T1 = 147	T0 = 27; T1 = 177
rcc/mOhm	0.6	0.7	0.7	0.8
rch/mOhm	0.6	0.7	0.7	0.8
Rth,c/K∙W^−1	31.0	31.5	30.9	31.0
Rth,h/K∙W^−1	23.5	23.3	22.8	23.0

Table 2.

Thermal resistances and contact electrical resistances under different temperatures.

3.2. Test setup

System for measuring output performance of thermoelectric modules was established, mainly including electric heating plate controlled by PID, adjustable load, circulatory cooling unit, thermal imaging device, temperature and voltage data acquisition units, etc., with its basic structure as shown in Figure 2. Electric heating plate is used as heat source, with temperature control precision of ±0.1 K and temperature ranging from room temperature to 773 K. Cooling unit, which consists of heat sink, water tank, flow meter and flow valve, etc., takes cold water as the coolant. Heat sink is made of red-copper and its temperature could be adjusted by controlling flowrate of cooling water. In addition, some thermal conductive filler is pasted on both sides of module to reduce thermal resistance between module, heat source and heat sink. Electrical current in the circuit is obtained via measuring voltage on both ends of sampling resistor (metal film precision resistor: 0.2 Ohm, precision of ±1%).

Voltage and temperature signal are acquired by 9207 and 9214 acquisition card of National Instruments (NI) Company, with precision of ±0.5%. Data to be acquired include as follows: (1) temperature of heat source and heat sink; (2) temperature of the coolant (water) inside heat sink and water tank; (3) voltage on adjustable load and sampling resistor. K-type thermocouples with diameter of 1 mm are inserted in heat source and heat sink to measure temperature values. Actually, even though electric heating plate is controlled by PID, heat source temperature still fluctuates during the change of load resistance. In order to eliminate impacts of such transient effect, data shall not be acquired until the heat source and heat sink temperatures are stable.

Figure 2.

Configuration of the output performance test system for the thermoelectric modules.

Test of energy efficiency is not undertaken due to its complexity, where the heat flow into the hot side of the module must be measured or evaluated. An effective way is to adopt heat flux sensor and bury it just under the module. But that would impact heat conduction between heat source and the module, leading to higher thermal resistance. And heat flux sensors of high temperature enduring are really costly. Another useful method is by calculating electrically generated heat in heat source, and at the same time, radiation and convective heat loss must be subtracted, as is introduced in Ref. [10].

3.3. Comparison of test results with calculation

Figures 3 and 4 show variations of output power P_out with load R_L at four temperature conditions, acquired by physical model calculation, ANSYS simulation and experiment. ANSYS method will be introduced in the next part. Figure 3 shows data at heat sink temperature T₀ = 300 K, while Figure 4—at T₀ = 296 K. Figures 5 and 6 are the corresponding current-voltage (I-V) characteristics. R_L results are disposed in the same way. From the results, it is found, that the output power has maximum value with the increase of load. And current is linearly related to voltage. Calculation results are well coincident with ANSYS results, and they are both a little higher than experimental data. Under the four temperature conditions, values of maximum output power are 2.5, 2.6, 2.8 and 1.1% higher than experimental results with T₁ changing from high to low. They are especially coincident well, when temperature difference ΔT is small. From the analysis follows, deviation of calculated results is caused mainly by taking thermal conductivity and electrical resistivity as constant (i.e., using the mean values), when solving Eq. (1). When ΔT is small, then material physical parameters vary within a narrow range. So, values of parameters are close to the real values. Otherwise, when ΔT is large, material physical parameters change within a large scale, leading to a great deviation of calculations.

Figure 3.

Dependences of output power on load resistance: calculations, experiments and ANSYS at T₀ = 300 K.

Figure 4.

Dependences of output power on load resistance: calculations, experiments and ANSYS at T₀ = 296 K.

Figure 5.

I-V characteristics of thermoelectric module: calculations, experiments and ANSYS at T₀ = 300 K.

Figure 6.

I-V characteristics of thermoelectric module: calculations, experiments and ANSYS at T₀ = 296 K.

4. Introduction to TEG simulation in ANSYS

4.1. TEG cell model

Figure 7.

The geometry of TEG cell in ANSYS and its mesh.

By software simulation, TEG performance can be achieved both in thermal and in electrical aspects. But it is not direct to cognize and understand the influence of thermoelectric effects, when compared with the above physical model. In this part, TEG cell model is set up by ANSYS, and geometry and meshing methods are illustrated in Figure 7. Thickness and cross-sectional area of thermoelements are 1.6 mm and 1.4 mm × 1.4 mm, respectively. Other geometry parameters are shown in Figure 7. Thermoelectric module consists mainly of p-n thermoelements, current-conducting copper straps and ceramic substrates for heat conducting and electric insulation. Thermoelements and copper strap are meshed by element SOLID226 in ANSYS. This type of element contains 20 nodes with voltage and temperature as the degrees of freedom. It can simulate 3D thermal-electrical coupling field. Element SOLID90 is used to mesh ceramic substrate. It has 20 nodes with temperature as the degree of freedom. Load resistance is simulated by element CIRCU124.

Contact properties of the leg-strap junction are implemented with element pairs CONTACT174/TARGET170. Detailed finite element formulations in ANSYS are introduced in [4], and the range of contact thermal conductivity and electrical resistivity is explicated in [11].

4.2. APDL codes for TEG simulation

ANSYS Parametric Design Language (APDL) is widely used for programed simulation. The following APDL codes have taken temperature variation of materials properties, thermal contact and thermal radiation (although its influence is very weak) into consideration. According to the practical requirements, the readers could use the code more concisely by neglecting certain physical effects. The unit referring to length is meter and the temperature unit is Celsius.

! defining the TEG cell dimensions

ln=1.6e-3 ! n-type thermoelement thickness

lp=1.6e-3 ! p-type thermoelement thickness

wn=1.4e-3 ! p-type thermoelement width

wp=1.4e-3 ! p-type thermoelement width

d=1.0e-3 ! Distance between the thermoelements

hs=0.2e-3 ! copper strap thickness

hc=1e-3 !substrate thickness

! definition of several physical parameters

rsvx=1.8e-8 ! copper electrical resistivity

kx=200 ! copper thermal conductivity

kxs=24 !substrate thermal conductivity

T1=250 ! temperature of heat source

T0=30 ! temperature of heat sink

Toffst=273 ! temperature offset

! defining TEG output parameters and the load

*dim,P0,array,1 ! defining P0 as the output power

*dim,R0,array,1 ! defining R0 as the load

*dim,Qh,array,1 ! defining Qh as the heat flow into the TEG cell

*dim,I,array,1 ! defining I as the current

*dim,enta,array,1 ! defining enta as the energy efficiency

*vfill,R0(1),ramp,0.025 ! setting the load (Ohm)

! pre-processing before calculation, defining element type, building the structure and meshing

/PREP7

toffst,Toffst ! set temperature offset

et,1,226,110 ! 20-node thermoelectric brick element

et,2,shell57 ! shell57 element for radiation simulation

et,3,conta174 ! conta174 element for contact simulation

et,4,targe170 ! target170 element for contact simulation

keyopt,3,1,4 ! taking temperature and voltage as the degree of freedom

keyopt,3,9,0

keyopt,3,10,1

keyopt,4,2,0

keyopt,4,3,0

! Temperature data points

mptemp,1,25,50,75,100,125,150

mptemp,7,175,200,225,250,275,300

mptemp,13,325,350

! Seebeck coefficient of the n-type material (V·K⁻¹)

mpdata,sbkx,1,1,-160e-6,-168e-6,-174e-6,-180e-6,-184e-6,-187e-6

mpdata,sbkx,1,7,-189e-6,-190e-6,-189e-6,-186.5e-6,-183e-6,-177e-6

mpdata,sbkx,1,13,-169e-6,-160e-6

! electrical resistivity of the n-type material (Ohm*m)

mpdata,rsvx,1,1,1.03e-5,1.06e-5,1.1e-5,1.15e-5,1.2e-5,1.28e-5

mpdata,rsvx,1,7,1.37e-5,1.49e-5,1.59e-5,1.67e-5,1.74e-5,1.78e-5

mpdata,rsvx,1,13,1.8e-5,1.78e-5

! thermal conductivity of the n-type material (m* K⁻¹)

mpdata,kxx,1,1,1.183,1.22,1.245,1.265,1.265,1.25

mpdata,kxx,1,7,1.22,1.19,1.16,1.14,1.115,1.09

mpdata,kxx,1,13,1.06,1.03

! Seebeck coefficient of the p-type material (V·K⁻¹)

mpdata,sbkx,2,1,200e-6,202e-6,208e-6,214e-6,220e-6,223e-6

mpdata,sbkx,2,7,218e-6,200e-6,180e-6,156e-6,140e-6,120e-6

mpdata,sbkx,2,13,101e-6,90e-6

! electrical resistivity of the p-type material (Ohm*m)

mpdata,rsvx,2,1,1.0e-5,1.08e-5,1.18e-5,1.35e-5,1.51e-5,1.7e-5

mpdata,rsvx,2,7,1.85e-5,1.98e-5,2.07e-5,2.143e-5,2.15e-5,2.1e-5

mpdata,rsvx,2,13,2.05e-5,2.0e-5

! thermal conductivity of the p-type material (m* K⁻¹)

mpdata,kxx,2,1,1.08,1.135,1.2,1.25,1.257,1.22

mpdata,kxx,2,7,1.116,1.135,1.13,1.09,1.12,1.25

mpdata,kxx,2,13,1.5,2.025

! material property for cooper strap

mp,rsvx,3,rsvx

mp,kxx,3,kx

! material property for the substrate

mp,kxx,4,kxs

!radiation property for the p-n materials

mp,emis,5

! contact friction coefficient

mp,mu,6,0

! build the TEG cell structure

block,d/2,wn+d/2,-ln,0,,t

block,-(wp+d/2),-d/2,-lp,0,,t

block,d/2,wn+d/2,,hs,,t

block,-(wp+d/2),-d/2,,hs,,t

block,-d/2,d/2,,hs,,t

block,-(wp+d/2),-d/2,-lp,-(lp+hs),,t

block,d/2,wn+d/2,-ln,-(ln+hs),,t

block,-(wp+d/2),wn+d/2,hs,hs+hc,,t

block,-(wp+d/2),wn+d/2,-(lp+hs),-(lp+hs+hc),,t

! glue the copper strap and the substrate

vsel,s,loc,y,0,hs

vsel,a,loc,y,hs,hc+hs

vglue,all

allsel

vsel,s,loc,y,-lp-hs,-lp

vsel,a,loc,y,-lp-hs-hc,-lp-hs

vglue,all

allsel

! meshing the TEG cell structure

numcmp,all

mshape,0,3d

mshkey,1

type,1

mat,3

lsel,s,loc,x,-d/2,d/2

lsel,r,loc,y,0

lsel,r,loc,z,t

lesize,all,d/3

vsel,s,loc,x,-d/2,d/2

vsel,r,loc,y,0,hs

vsweep,all

allsel

esize,ww/3

type,1

mat,3

vsel,s,loc,y,0,hs

vsel,u,loc,x,-d/2,d/2

vsweep,all

vsel,s,loc,y,-lp-hs,-lp

vsweep,all

type,1

mat,1

vsel,s,loc,x,d/2,d/2+wn

vsel,r,loc,y,-ln,0

vmesh,all

mat,2

vsel,s,loc,x,-(wp+d/2),-d/2

vsel,r,loc,y,-lp,0

vmesh,all

type,1

mat,4

vsel,s,loc,y,hs,hs+hc

vsel,a,loc,y,-lp-hs-hc,-lp-hs

vsweep,all

allsel

! defining the contact parameters

r,5 ! selecting the thermal contact conductivity and resistivity

RMORE,

rmore,,7e5 ! setting the thermal contact conductivity

rmore,0.67e8,0.5 ! setting the thermal contact resistivity

! defining the contact layer between p-leg and upper copper strap

vsel,s,loc,y,0,hs

asel,s,ext

asel,r,loc,y,0

nsla,s,1

nsel,r,loc,x,-(wp+d/2),-d/2

type,3

mat,6

real,5

esurf

allsel

! defining the target layer between p-leg and upper copper strap

vsel,s,mat,,2

asel,s,ext

asel,r,loc,y,0

nsla,s,1

type,4

mat,6

esurf

allsel

! defining the contact layer between n-leg and upper copper strap

vsel,s,loc,y,0,hs

asel,s,ext

asel,r,loc,y,0

nsla,s,1

nsel,r,loc,x,d/2,d/2+wn

type,3

mat,6

real,5

esurf

allsel

! defining the target layer between n-leg and upper copper strap

vsel,s,mat,,1

asel,s,ext

asel,r,loc,y,0

nsla,s,1

type,4

mat,6

esurf

allsel

! defining the contact layer between p-leg and bottom copper strap

vsel,s,loc,y,-hs-lp,-lp

vsel,r,loc,x,-wp-d/2,-d/2

asel,s,ext

asel,r,loc,y,-lp

nsla,s,1

type,3

mat,6

real,5

esurf

allsel

! defining the target layer between p-leg and bottom copper strap

vsel,s,mat,,2

asel,s,ext

asel,r,loc,y,-lp

nsla,s,1

type,4

mat,6

esurf

allsel

! defining the contact layer between n-leg and bottom copper strap

vsel,s,loc,y,-hs-ln,-ln

vsel,r,loc,x,d/2,d/2+wn

asel,s,ext

asel,r,loc,y,-ln

nsla,s,1

type,3

mat,6

real,5

esurf

allsel

! defining the target layer between n-leg and bottom copper strap

vsel,s,mat,,1

asel,s,ext

asel,r,loc,y,-ln

nsla,s,1

type,4

mat,6

esurf

allsel

! defining the shell element for radiation simulation, outputting radiation matrix

! defining the shell element for copper strap

type,2

aatt,3,,2

asel,s,loc,x,-(wp+d/2),wn+d/2

asel,r,loc,y,0,hs

asel,u,loc,y,0

asel,u,loc,y,hs

amesh,all

allsel

asel,s,loc,x,-d/2,d/2

asel,r,loc,y,0

amesh,all

allsel

aatt,3,,2

asel,s,loc,x,-(wp+d/2),wn+d/2

asel,r,loc,y,-lp-hs,-lp

asel,u,loc,y,-lp

asel,u,loc,y,-lp-hs

amesh,all

allsel

aatt,4,,2

asel,s,loc,x,-d/2,d/2

asel,r,loc,y,-lp-hs

amesh,all

! defining the shell element for p-n thermoelements

allsel

aatt,5,,2

asel,s,loc,x,-(wp+d/2),wn+d/2

asel,r,loc,y,-lp,0

asel,u,loc,y,-lp

asel,u,loc,y,0

amesh,all

! defining the space node for radiation simulation

n,10000,0,0,3e-3

fini

! using radiation matrix method

/aux12

emis,3,1 ! setting the emissivity

emis,4,1

emis,5,1

allsel

geom,0

stef,5.68e-8 ! setting the Stefan-Boltzmann constant

vtype,hidden

space,10000

write,teg,sub ! outputting the radiation super element

fini

/prep7

! deleting the shell elements and the corresponding mesh

allsel

asel,s,type,,2

aclear,al

etdele,2

allsel

et,5,matrix50,1 ! defining radiation matrix element

! defining boundary conditions and the load

nsel,s,loc,y,hs+hc ! TEG cell hot side

cp,1,temp,all ! coupling of temperature degree of freedom

nh=ndnext(0) ! getting the master node

d,nh,temp,Th ! setting the temperature constraint to the hot side

nsel,all

nsel,s,loc,y,-(ln+hs+hc) ! selecting the TEG cell cold side

d,all,temp,Tc ! setting the temperature constraint to the cold side

nsel,s,loc,y,-(ln+hs),-ln

nsel,r,loc,x,d/2+wn

cp,3,volt,all ! electrical coupling

nn=ndnext(0) ! getting the master node

d,nn,volt,0 ! setting the ground connection node

nsel,all

nsel,s,loc,y,-(lp+hs),-lp

nsel,r,loc,x,-(wp+d/2)

cp,4,volt,all ! ! electrical coupling

np=ndnext(0) ! getting the master node

nsel,all

type,5

allsel

d,10000,temp,300 ! setting the temperature of the space node

se,teg,sub ! reading the radiation super element

et,6,CIRCU124,0 ! setting the load resistor element

fini

/prep7

! setting the load value and property

r,1,R0(1)

type,6

real,1

numcmp,all

e,np,nn

esel,s,type,,6

circu_num=elnext(0) !getting circuit element number

allsel

fini

! starting the calculation

/SOLU

antype,static ! solution type

cnvtol,heat,1,1.e-3 ! setting the converging value for heat condition

cnvtol,amps,1,1.e-3 ! setting the converging value for the current

neqit,50 ! calculation iteration step

solve ! starting solving

fini

*get,P0(1),elem,circu_num,nmisc,1 ! getting the output power of the TEG cell

*get,Qh(1),node,nh,rf,heat ! getting the heat flow into the TEG cell

*get,I(1),elem,circu_num,smisc,2 ! getting the current

*voper,enta,P0,div,Qh ! calculating the energy efficiency of the TEG cell

5. Conclusions

The built one-dimensional model, which is validated by test results, can calculate TEG output power and energy efficiency accurately. By simplifying this model, it is convenient to analyze influences of different thermal and electrical parameters on TEG performance. And basic factors to enhance TEG output power and energy efficiency are extracted. At last, ANSYS simulation considering thermal contact and radiation effects for TEGs is introduced briefly, and basic APDL codes are shared.