The TEM model is illustrated in Figure 1. The hot/cold substrates and copper pads are required to manufacture the TEM, yet only the TE legs are considered in the present work. This implies that thermal resistances caused by substrates and copper pads, thermal/electrical contacts, and convective/radiative thermal leaks inside the TEM are not considered. The abovementioned source terms, which are thermal resistances/contacts and various thermal leaks, can be conveniently incorporated. The present work focuses on solving the highly non-linear energy equation when considering temperature-dependent material properties. Note that the generated current results in further heat pumping from the hot-end to the cold-end when the TEM is subjected to the power generation mode. This is confirmed in various experiments (Li et al., 2019). The z-coordinate in Figure 1B starts from the cold-end of the TE leg and directs to the hot-end.

The governing energy equation for the TE leg can be derived on the basis of energy conservation, the first Thomson relation and second Thomson relation. The governing energy equation can be expressed as follows after treating the current (I) as scalar. d − k T ⋅ d T d z d z = I 2 ρ A 2 ± I A T d α d T d T d z , (1) where + is for P-type leg and − is for N-type leg. Note that this governing equation is consistent with previous references (Ju et al., 2017; Lee et al., 2018). This partial equation can be solved under the assumption that the properties (α, ρ, and k) of the TE material are independent of the temperature. The obtained result is as follows: T = − I 2 ρ 2 k A 2 z 2 + T h − T c L + I 2 ρ L 2 k A 2 z + T c , (2) where T _h, T _c, L, and A are the hot-end temperature, cold-end temperature, length of the TE leg, and cross-sectional area of the TE leg, respectively. The temperature distribution along the TE leg is parabolic on the basis of Eq. 2. Therefore, Eq. 2 is referred to as the second-order temperature distribution. In case that temperature-dependent material properties are considered, the temperature distribution should be more complex.

Bi₂Te₃-based TE material is used widely and can be obtained easily in the open market. Therefore, temperature-dependent material properties of the Bi₂Te₃-based TE material from Sagreon Co., Ltd., China, are used in the present work. However, the proposed model is not limited to the Bi₂Te₃-based TE material. Third-order polynomial fittings are applied (Li et al., 2020b) as follows: k p T = t p 1 + t p 2 T + t p 3 T 2 + t p 4 T 3 , (3-1) k n T = t n 1 + t n 2 T + t n 3 T 2 + t n 4 T 3 , (3-2) ρ p T = e p 1 + e p 2 T + e p 3 T 2 + e p 4 T 3 , (4-1) ρ n T = e n 1 + e n 2 T + e n 3 T 2 + e n 4 T 3 , (4-2) α p T = s p 1 + s p 2 T + s p 3 T 2 + s p 4 T 3 , (5-1) α n T = s n 1 + s n 2 T + s n 3 T 2 + s n 4 T 3 . (5-2) The constants in Eqs 3–5 are as follows: t p 1 = 4.389 , t p 2 = − 1.8168 ✕ 10 − 2 , t p 3 = 2.437956 ✕ 10 − 5 , t p 4 = 4.793196 ✕ 10 − 10 , t n 1 = 4.09878 , t n 2 = − 1.4976 ✕ 10 − 2 , t n 3 = 1.799196 ✕ 10 − 5 , t n 4 = 1.692996 ✕ 10 − 9 , e p 1 = − 6.7074 ✕ 10 − 6 , e p 2 = 5.09 ✕ 10 − 8 , e p 3 = 6.33243 ✕ 10 − 11 , e p 4 = − 5.31761 ✕ 10 − 14 , e n 1 = − 1.51744 ✕ 10 − 5 , e n 2 = 1.142 ✕ 10 − 7 , e n 3 = − 8.17056 ✕ 10 − 11 , e n 4 = − 5.18487 ✕ 10 − 15 , s p 1 = − 1.0915819 ✕ 10 − 4 , s p 2 = 1.67585 ✕ 10 − 6 , s p 3 = − 2.12 ✕ 10 − 9 , s p 4 = 4.43743 ✕ 10 − 14 , s n 1 = − 4.3833365 ✕ 10 − 4 , s n 2 = 2.90422 ✕ 10 − 6 , s n 3 = − 9.76 ✕ 10 − 9 , s n 4 = 1.01202 ✕ 10 − 11 .

For a control volume (CVM) of AΔz (Δz is the length of the control volume), number the control volume from the cold-end to the hot-end, and the first control volume has two nodes, i.e., i = 0 and i = 1. An algorithm to calculate the temperature of the (i+1)^th control volume can be derived on the basis of energy conservation and Eq. 1 as follows: k T i A T i − T i − 1 Δ z = k T i A T i + 1 − T i Δ z + I 2 ρ T i Δ z A ± μ T i I A T i + 1 − T i − 1 2 Δ z A Δ z , (6)

Note that the (i + 1)^th control volume is physically realistic. Such numbering rules avoid the confusing zero control volume, and it is beneficial during coding. The second and third terms in the right-hand side of Eq. 6 represent the Joule heat and Thomson heat, respectively. + is for P-type TE leg and − is for N-type TE leg. As a result, Eq. 6 can be combined to obtain the temperature distribution along the TE leg based on the CVM as follows: T i = 0.5 T i − 1 + T i + 1 + I 2 ρ T i Δ z 2 2 k T i A 2 ± μ T i I Δ z 4 k T i A T i + 1 − T i − 1 (7)

On the other hand, by applying a second-order central finite difference algorithm on Eq. 1, the following algorithm can be derived to obtain the temperature distribution based on the finite difference method (FDM) as follows: T i = 0.5 T i − 1 + T i + 1 + I 2 ρ T i Δ z 2 2 k T i A 2 ± μ T i I Δ z 4 k T i A T i + 1 − T i − 1 + k d T i 8 k T i T 2 i + 1 + T 2 i − 1 − 2 T i + 1 T i − 1 (8) where + is for P-type TE leg and − is for N-type TE leg. k _d is the first-order derivative of k. By comparing Eqs 7 and 8, these two numerical models were found to be the same, except for the fourth term in the right-hand side of Eq. 8. This term represents the influence of the local change rate of thermal conductivity on temperature distribution, whereas the thermal conductivity inside a CVM is assumed to be unchanged.

The electric current and internal electrical resistances are coupled with the temperature distribution. As a consequence, Eqs 7 and 8 become difficult to be solved. An iteration method is required to decouple the abovementioned problem. The open circuit voltage can be obtained as follows: V O C = N ∫ T c T h α P − α N d T = N s p n 1 T h − T c + 1 2 s p n 2 T h 2 − T c 2 + 1 3 s p n 3 T h 3 − T c 3 + 1 4 s p n 4 T h 4 − T c 4 ] , (9) where spn ₁ = sp ₁ − sn ₁, spn ₂ = sp ₂ − sn ₂, spn ₃ = sp ₃ − sn ₃, and spn ₄ = sp ₄ − sn ₄. The internal electrical resistance cannot be obtained because the temperature distribution is not known. However, the averaged internal electrical resistance can be used as an initial value to decouple Eqs 7 and 8 as follows: R i n = N A ∫ 0 L ρ p + ρ n d z = N A ρ P a v + ρ N a v L , (10) where N is the number of TE couples inside the TEM. The average electrical resistivity can be determined based on Eq. 4 as follows: ρ P a v = 1 T h − T c ∫ T c T h ρ p T d T = 1 T h − T c e p 1 T h − T c + 1 2 e p 2 T h 2 − T c 2 + 1 3 e p 3 T h 3 − T c 3 + 1 4 e p 4 T h 4 − T c 4 ] , (11-1) ρ N a v = 1 T h − T c ∫ T c T h ρ n T d T = 1 T h − T c e n 1 T h − T c + 1 2 e n 2 T h 2 − T c 2 + 1 3 e n 3 T h 3 − T c 3 + 1 4 e n 4 T h 4 − T c 4 ] , (11-2)

Hence, the electric current can be determined when the external load (R _ex) is provided as follows: I = V O C R i n av + R ex , (12)

The TE efficiency is obtained when R _ex = R _in and is defined as follows: ζ = P max Q h = V O C 2 Q h R i n , (13) where Q _h is the heat flow from the hot-end of the TEM, and it is calculated as follows: Q h = N − k P T h , P A d T P d z z = L + α P T h , P I T h , P − k N T h , N A d T N d z z = L + α N T h , N I T h , N ] (14)

Boundary conditions, which are T = T _c at z = 0 m and T = T _h at z = L m have to be employed in the numerical model. Note that only constant temperature conditions are employed in the present study because most of the previous studies on TEG have applied this type of boundary condition. For the constant heat flux boundary condition, a recent insightful investigation has reported observations in detail (Min, 2022). Eqs 2–14 were coded with C language to explore the influence of temperature-dependent material properties on TEG performance and to find out the underlying thermal–electrical multi-physics during TE conversion.

As stated previously, the averaged internal electrical resistance is used as an initial value to decouple Eqs 7 and 8. Then, the new internal electrical resistance must be re-calculated once the temperature distribution is obtained. Hereby, iterations have to be performed to obtain new temperature distributions until the internal electrical resistance between two iterations becomes small enough (e.g., 0.0001 Ω). In general, several iterations (less than five iterations) are enough. This also implies that a satisfied temperature distribution is obtained. The CPU time for each set of input parameters (a, L, F, T _h, T _c, and R _ex) is 3.82 s with a 12th Intel i7-1260P processor, and the total CPU time for structure optimization depends on the number of parameter sets. Figure 2 shows the flow chart of the developed numerical model. Note that another subroutine is required to obtain the leg temperature distribution, T _P (z) and T _N (z) in subroutine B. Hence, there are three different levels of the loop program. An ordinary least squares method was used to obtain T _P (z) and T _N (z), and this subroutine is not shown for the sake of briefness.

Figure 3 shows the experimental setup for the verification of the developed numerical model. The tested TEM is a commercially available TEG1-12708 from Sagreon Co., Ltd., China. The dimensions are 40 mm × 40 mm × 3.8 mm, and the temperature-dependent properties can be found in our previous work (Li et al., 2020b). In order to obtain a constant hot-end temperature, high temperature silicon oil with a boiling temperature of 573 K was used to transfer heat from the electrical heater to the copper-based heat collector and distributor. The hot-end temperature was set with a temperature controller. The cold-end temperature was ensured with an alumina heat sink, a blower, and a speed controller. Thermal grease and a weight of 10 kg were applied to decrease the thermal contacts between the TEM and heat-distributor/heat-sink. The generated power was measured with an electronic load (Prodigit 3311F) by adjusting different external loads. The accuracy of the electronic load was ±0.5%. Two thermocouples with diameters of 0.5 mm were carefully installed in the pre-grooved heat-distributor/heat-sink. The accuracy of the thermocouple was ±0.5%. All experimental cases were run at least 2 h before performing measurements so as to eliminate possible unsteady conditions. The cold-end temperature was maintained at 323 K by adjusting the blower speed. The cold-end temperature of 323 K was necessary for forced air cooling because a certain temperature difference between the cold-end temperature and atmospheric temperature is required. The atmospheric temperature during the experiments was 299 K.

As shown in Figure 4, the predicted results agree well with the experimental data under the same hot-/cold-end temperatures. The average experimental electric powers are 0.023 W, 0.107 W, 0.238 W, and 0.412 W lower than that obtained by the developed numerical model under the hot-end temperatures of 373 K, 423 K, 473 K, and 523 K, respectively. The reasons for the error include the negligence of thermal/electrical contacts and the one-dimensional assumption of the TE leg. For example, the effective temperature difference is lower than that in the numerical model, which leads to less electric power than that by the numerical model. Hence, the developed numerical model can predict the TEG performance with acceptable errors. Furthermore, the developed numerical model can be used to reveal various aspects of the underlying thermal–electrical multi-physics as given in the following sections.

Figure 5 shows the comparisons on the convergent performance between the CVM and FDM in obtaining T _P (z) and T _N (z) in subroutine B. In order to present the convergent process, T _p − T _2nd is used as the monitor parameter. As shown in Figure 5, the FDM performs much better than the CVM, that is, the FDM is found to have converged much faster with fewer grids than the CVM. For the FDM, it usually takes several seconds to reach satisfied solutions of T _P (z) and T _N (z). Note that the maximum T _p − T _2nd is found to be approximately 0.5 K when the CVM is employed, yet the corresponding value is approximately 1.3 K for the FDM. The underlying reason for the abovementioned phenomenon is the inclusion of the local change rate of thermal conductivity in the FDM, that is, the fourth term in the right-hand side of Eq. 8. This reveals that the local change rate of thermal conductivity considerably affects the temperature distribution. In fact, the thermal conductivity is of more priority than the Seebeck coefficient and electrical resistivity during the optimization of TEG performance (Shen et al., 2021). Hence, the local change rate of thermal conductivity inevitably affects the temperature distribution.

Figure 6 presents the temperature difference for P-/N-type legs with different numerical methods. As shown in Figure 6, the maximum temperature differences between the FDM and CVM can be as large as 1.0 K, but the locations for the abovementioned maximum temperature difference are different. For the P-type leg, the abovementioned location is near the hot-end of TEM, whereas it is near the cold-end of TEM for the N-type leg. To the best of the authors’ knowledge, this phenomenon has not been reported in previous studies, yet our experiments cannot confirm this finding due to the limitation of measurement techniques. Further theoretical studies and advanced experiments have to be carried out to verify this finding. Therefore, the temperature distributions for the P-/N-type legs are found to be different, which are presented and discussed in detail in the following section. In the following sections, all results are obtained on the basis of the FDM.

Figure 7 presents the temperature distribution of the P-/N-type legs when compared to the linear and parabolic temperature distributions. A well-accepted conclusion for the second-order temperature distribution is that the location of the maximum temperature difference is at the central point of the leg (Rowe and Gao, 1998). However, temperature distributions considering temperature-dependent material properties are inconsistent with the second-order temperature distribution. The temperature of the P-type leg is always larger than that predicted by the second-order temperature distribution in every location of the leg, whereas the temperature of the N-type leg is always lower than that predicted by the second-order temperature distribution. The maximum temperature difference for the P-type leg is approximately 2.2 K when compared to that of linear temperature distribution, whereas the corresponding value is approximately 0.9 K for the second-order temperature distribution. Moreover, the location of the maximum temperature difference for the P-type leg is not at the central point of the leg. The abovementioned finding answers the doubt of the influence degree of temperature-dependent material properties on the temperature distribution. It is found that temperature-dependent material properties cause an unexpected influence on temperature distributions, and the influence degree is more serious than that between the second-order and linear temperature distributions.

The developed numerical model can predict the response of temperature distribution under different external loads, which is shown in Figure 8. It is found that the ratio of the external load to internal electrical resistance (R _ex/R _in) considerably affects the temperature distribution. A conclusion can be obtained from Figure 8 by revising Eqs 1 and 12: the electric current notably affects the temperature distribution. A low R _ex/R _in results in a large current. Therefore, augmented Peltier heat pumping occurs, which leads to different temperature distributions. Note that the Thomson effect is also involved, and it is difficult to distinguish the influence level by the Thomson effect in Figure 8. However, more detailed results and discussions are presented in the following section related to this concern.

An interesting exploration is to intentionally remove the Thomson term from Eq. 8, and the results are shown in Figure 9 and Figure 10. As shown in Figure 9, the temperature difference caused by the Thomson effect for the P-/N-type legs are considerable, and R _ex/R _in also plays an important role in causing the abovementioned temperature difference. By comparing Figures 6, 7, and 9, it is found that the influence of the local change rate of thermal conductivity and the Thomson effect in the P-type leg are in the same phase, which enlarges the temperature difference in Figure 7. However, the influence of the local change rate of thermal conductivity is out of phase with the influence of the Thomson effect in the N-type leg. Hence, the temperature difference of the N-type leg is smaller than that of the P-type leg as seen in Figure 7.

The heat absorption and release by the Thomson effect are quantitatively shown in Figure 10. It is clearly seen that both heat absorption and release occur in the P-/N-type legs. That is, heat is released in a certain length of the leg near the cold-end of the TEM, whereas heat is absorbed in the remaining length of the leg near the hot-end of the TEM. The abovementioned phenomenon reveals the complicated nature of TE conversion.

The electric power (P) and TE efficiency (ξ) under different hot-end temperatures and R _ex/R _in are shown in Figure 11. Several sets of data are extracted from Figure 11A to compare with the experimental data as shown in Figure 4. For a constant cold-end temperature (323 K), the electric power is mainly determined by the hot-end temperature, and R _ex/R _in is another important parameter. A new finding shown in Figure 11A is that the obtained curve surface bends with a larger amplitude as the hot-end temperature increases. On the other hand, the TE efficiency monotonously increases with the temperature difference when constant material properties are assumed. However, the abovementioned monotonous increase is no longer existent when temperature-dependent material properties are considered. As shown in Figure 11B, the maximum TE efficiency is found to be near the hot-end temperature of 500 K and further increasing the hot-end temperature leads to lower TE efficiency than that near 500 K even though the electric power continues to increase, as shown in Figure 11A. The abovementioned phenomenon is caused by the fact that the thermal conductivity of the TE material monotonously increases with temperature, but the Seebeck coefficient exhibits a parabolic shape with temperature. Note that the TE efficiency could be over predicted because thermal resistances caused by substrates and copper pads, thermal/electrical contacts, and convective/radiative thermal leaks inside the TEM are not considered, which is stated in Section 2.1. However, this will not alter the conclusion of the present work.

Figure 12 shows the structural optimization results concerning electric power and TE efficiency. The electric power can be customized according to particular applications because the electric power is sensitive to all structural parameters (a, L, and F). Note that the results in Figure 12 are obtained under the assumption that the thermal leaks by air conduction and thermal radiation inside the TEM are neglected. Minor influence on electric power is expected because the ratio of the thermal leak to the total heat flux passing through the TEM is only several percentages even though changing F results in different quantities of thermal leaks. Hence, Figure 12 presents valuable optimization results in the aspect of structure configuration. The ordinary method of predicting the electric power is to assume the temperature-dependent material properties as constants, i.e., averaged material properties between hot-/cold-ends. This assumption is effective, and only minor differences can be found when using the present numerical model, which is shown in Figure 12B. It is found that the electric power differences (P _max − P _maxav) are within 1% of the developed numerical model and the ordinary model with averaged material properties. P _maxav is the maximum electric power with average TE material properties.

TE efficiency remains unchanged under different structural parameters, as shown in Figure 12C. This reveals that an attempt to augment the TE efficiency with the same TE material through structural optimization has little potential. Note that the TE efficiency was found to be changed slightly when modifying the TEM structure (Zhao et al., 2021). This phenomenon is caused by the thermal leak inside the TEM (Lee et al., 2018). In case that Eq. 14 includes the conduction heat loss and radiative thermal loss between the hot-/cold-ends (Lee et al., 2018), the TE efficiency difference (ξ − ξ _lk) versus structural parameters is shown in Figure 12D, where ξ _lk is the TE efficiency when thermal leaks inside the TEM are considered in Eqs 13 and 14. As shown in Figure 12D, the major parameter affecting the TE efficiency is the occupancy ratio, and the leg length also plays a moderate influence on TE efficiency. The influence level by thermal leaks could reach up to 15%. A low occupancy ratio leads to a serious downgrade of TE efficiency. The influence of thermal leaks on the TE efficiency is not the major task of the present work.