Analyzing the Performance of Thermoelectric Generators with Inhomogeneous Legs: Coupled Material–Device Modelling for Mg₂X-Based TEG Prototypes

Abstract

Thermoelectric generators (TEGs) possess the ability to generate electrical power from heat. As TEGs are operated under a thermal gradient, inhomogeneous material properties—either by design or due to inhomogeneous material degradation under thermal load—are commonly found. However, this cannot be addressed using standard approaches for performance analysis of TEGs in which spatially homogeneous materials are assumed. Therefore, an innovative method of analysis, which can incorporate inhomogeneous material properties, is presented in this study. This is crucial to understand the measured performance parameters of TEGs and, from this, develop means to improve their longevity. The analysis combines experimental profiling of inhomogeneous material properties, modelling of the material properties using a single parabolic band model, and calculation of device properties using the established Constant Property Model. We compare modeling results assuming homogeneous and inhomogeneous properties to the measurement results of an Mg₂(Si,Sn)-based TEG prototype. We find that relevant discrepancies lie in the effective temperature difference across the TE leg, which decreases by ~10%, and in the difference between measured and calculated heat flow, which increases from 2–15% to 9–16% when considering the inhomogeneous material. The approach confirms additional resistances in the TEG as the main performance loss mechanism and allows the accurate calculation of the impact of different improvements on the TEG’s performance.

1. Introduction

Thermoelectric (TE) materials can convert heat flow into electrical power. They have gained a lot of interest over the past decades as a green source of electrical energy [1] since 60% of the primary energy is lost as waste heat [2]. TE devices have the advantages of being lossless scalable, and negligible maintenance requirements due to their lack of moving parts. This has made them a reliable energy source in demanding fields, such as the aerospace industry, in which RTGs (Radioisotope Thermoelectric Generators) were used in several space missions [3]. Terrestrial applications are also being considered and researched, such as in the automotive industry and industrial processes [3,4,5], wearable medical devices [6,7,8], mobile storage of pharmaceuticals, and electronic devices [3].

A thermoelectric generator (TEG) is a device in which n- and p-type TE legs are connected electrically in series and thermally in parallel to generate an electrical current from a heat flow. Over the last decades, the optimization of the TE properties of various material classes has been the main focus of research as the first, very challenging step in the development chain of a TEG [9,10,11,12,13,14,15,16,17,18]. As a consequence, research on contacting techniques and on TE module building remained relatively scarce, and only TEGs based on (Bi,Sb)₂(Te,Se)₃ and PbTe/TAGS have reached commercial maturity. For many other promising material classes such as half-Heusler materials, Skutterudites, PbTe, and Zintl- phases, prototypes have been demonstrated [19,20,21,22,23]; partially, with limited stability due to an early development stage.

Another promising material class is solid solutions from Mg₂Si-Mg₂Sn, due to the high performance of the n-type material, the recent improvement of the p-type material, and the abundance and non-toxicity of its components [12,13,24,25,26]. It is also lightweight; which, is advantageous for mobile applications. For this material class, few prototypes have been reported recently [27,28], and we have performed the first efficiency measurement on an Mg₂(Si,Sn)-based TEG, demonstrating 3.6% conversion efficiency for a hot side temperature of 400 °C and a cold side temperature of 25 °C [29].

Locally resolved measurements showed a spatially varying Seebeck coefficient, indicating an inhomogeneous carrier concentration in the initially homogeneous material. As the device properties depend on the figure of merit $zT=\frac{{\alpha }^2\sigma }{\kappa }T$ of the employed material and the material properties α (Seebeck coefficient), σ (electrical conductivity), and κ (thermal conductivity), all depend on the carrier concentration that has a direct impact on the device’s performance. In fact, spatially inhomogeneous materials are quite a common feature in TEGs. With respect to our prototype, Mg₂(Si,Sn) was previously shown to be sensitive to Mg evaporation at expected working temperatures [30,31] and in addition to interdiffusion with foreign elements, such as the metallization layers which are included in a TEG design [32,33]. Both mechanisms usually lead to spatially inhomogeneous changes in the carrier concentration due to changes in the intrinsic or extrinsic defect densities. Other material systems show similar chemical or thermal instabilities, leading to inhomogeneity. Such as: Mg₃Sb₂, which is also sensitive to Mg loss [34]; CoSb₃ and half-Heusler compounds, which are sensitive to Sb loss [35]; as well as (Bi,Sb)₂(Te,Se)₃, which is sensitive to corrosion [36]. Generally, due to high temperatures (processing, application) and the proximity of multiple chemical elements, doping defects can form in the TE materials, which alter their properties [32]. On the other hand, local variation of properties can also be intended and designed, as thermoelectric material properties are generally quite strongly temperature-dependent; in addition, strategies such as segmentation [25,37] or grading are employed for performance optimization [38].

Those local changes in TE properties can be challenging when it comes to TEG modeling, used to predict the performance of a TEG, and in providing a reference to evaluate and understand measured data. The Constant Property Model (CPM) is generally used in TEG calculations for its simplicity and ease [25,39,40,41]. For real materials, averages corresponding to the relevant temperature range/profile are employed, resulting in a relatively accurate prediction of TEG performance [32,42]. However, these averages are obtained from temperature-dependent data; which, are usually obtained by integral bulk measurements on material samples before device manufacturing, and are therefore only available for homogeneous materials. On the other hand, locally resolved properties of inhomogeneous samples are measured typically only at room temperature [43,44]. To be able to model inhomogeneous materials in general, but also to calculate correct average values for the CPM, there is a need for a material model that is both spatially and temperature dependent. The performance of inhomogeneous thermoelectric materials was already modeled using the Effective Media Theory [45,46,47]; however, those works are generally applied to (intended) nanostructured thermoelectric materials, i.e., typically a composite structure made of particles of a first material embedded in a matrix made of a second material. This is a different case than that addressed here, as these materials are homogeneous on a mesoscale/macroscale.

In this study, an innovative analytical method is developed combining experimental carrier concentration profiling and the Single Parabolic Band (SPB) model combined with the CPM/continuum theory: the CPIM (Constant Property (model) for Inhomogeneous Materials). This method is able to capture the inhomogeneity of thermoelectric materials on a macroscopic/mesoscopic scale and is employed exemplarily to the measurement results of an Mg₂(Si,Sn)-based TEG. The predictions considering the actual inhomogeneous properties of the n-type legs are compared to those assuming a homogeneous material and the experimental data. When considering the inhomogeneous material, the deviation between calculations and measurement of the heat flow increases above the measurement uncertainty, and higher thermal losses through the TEG/heat exchanger coupling are observed. However, for both conditions, the high difference observed between calculated and measured electrical resistance likely indicates crack formation; which, is a commonly observed degradation mechanism in TEGs.

2. Methods

Experimental:

This work deals with the same TEG as reported in [29], where the building process and characterization methods are described. In this study, the solid solution Mg₂Si_0.3Sn_0.7 was chosen as the chemical composition for both p- and n-type materials of the TEG. Both materials have similar mechanical properties [48] and have a reported figure of merit among the best in the mid- to high-temperature range. Pictures and dimensions are reported in [29], while a schematic of the prototype is shown in Figure 1. The powder used to sinter the n-type material was synthesized and compacted, such as reported in [24]; the resulting properties are shown in Figure 2. The p-type powder is synthesized similarly, as described in [32].

The Seebeck coefficient and electrical conductivity were measured using an in-house device with a four-probe technique (HTSσ) [44,49]. The thermal conductivity $\mathrm{\kappa }$ was calculated from $\kappa =D\rho C_p$, where $D$, $\rho$, and $C_p$ are the sample thermal diffusivity, density, and heat capacity depending on the composition at constant pressure, respectively. $D$ was measured by a laser flash technique with a NETZSCH LFA 427 apparatus or with an XFA 467HT Hyperflash apparatus; ρ was measured using Archimedes’ method. $C_p=C_v^{\mathrm{D}\mathrm{P}}+\frac{9E_T^{2}T}{{\beta }_T\rho }$, where $C_v^{\mathrm{D}\mathrm{P}}$ is the Dulong-Petit limit; $E_T$ and ${\beta }_T$, respectively, are the linear coefficient of thermal expansion and isothermal compressibility dependent on composition [50]. The values for Mg₂Si_0.3Sn_0.7 are 2∙10⁻⁵ K⁻¹ and 2.07∙10⁻¹¹ Pa⁻¹, respectively [43]. The measurement uncertainties for $\alpha$, $\sigma$, and $\kappa$ are ±5%, ±5%, and ±8%, respectively, based on a comparison with the NIST low-temperature standard for the Seebeck coefficient [51] and internal reference measurements on a high-temperature standard [52]. Our estimates are comparable to those obtained in an international Round-robin test.

The metallization of the pellets was done using Al foils with Zn coatings, such as those reported in [53]. The resulting p-type legs are homogeneous while the n-type legs show a property gradient, presumably due to Zn diffusion.

CPIM (Constant Property model for Inhomogeneous Materials):

The CPIM relies on the application of the CPM (Constant Property Model) to inhomogeneous materials. The inhomogeneity is implemented in the model by experimentally obtaining a spatial distribution of the carrier concentration in the TE material, using a Seebeck coefficient microprobe and the SPB (Single Parabolic Band) model to link both quantities. The SPB is then used to obtain a spatial temperature-dependent distribution of all relevant quantities for the calculation of the average values used in the CPM. Each component of the model is adequately detailed below.

SPB model:

To capture the inhomogeneous properties of the n-type legs and be able to predict the properties at temperatures higher than room temperature, a single parabolic band model is employed. An SPB model allows us to calculate the macroscopic transport n-type properties based on a few material parameters: the reduced chemical potential ($\eta =\frac{E_F}{k_{\mathrm{B}}T}$), where $E_F$ is the Fermi energy and $k_{\mathrm{B}}$ is Boltzmann’s constant); the mobility parameters for acoustic phonon scattering (AP) and alloy scattering (AS) (${\mu }_{0,\mathrm{A}\mathrm{P}}$ and ${\mu }_{0,\mathrm{A}\mathrm{S}}$, respectively); and the density of states effective mass ($m_D^{\mathrm{*}}$) [50,54]. In this work, the lattice thermal conductivity (${\kappa }_{\mathrm{l}\mathrm{a}\mathrm{t}}$) is also used as an input parameter. The transport properties are obtained using the following equations; which, are given here in the specific case corresponding to AP and AS as relevant scattering mechanisms.

$$\alpha =\frac{k_{\mathrm{B}}}{e}\left(\frac{2F_1}{F_0}-\eta \right)$$

(1)

$$n=4\pi {\left(\frac{2m_D^{*}{k}_{\mathrm{B}}T}{h^2}\right)}^{1.5}{F}_{\frac{1}{2}}\left(\eta \right)$$

(2)

$${\mu }_{\mathrm{A}\mathrm{P}}={\mu }_{0,\mathrm{A}\mathrm{P}}·\psi \left(\eta \right)=\frac{\sqrt{8\pi }e{\hslash }^4\rho {v_l}^2}{3{E_{\mathrm{D}\mathrm{e}\mathrm{f}}}^2{m_s}^{2.5}{\left(k_{\mathrm{B}}T\right)}^{1.5}}\psi (\eta )$$

(3)

$${\mu }_{\mathrm{A}\mathrm{S}}={\mu }_{0,\mathrm{A}\mathrm{S}}·\psi \left(\eta \right)=\frac{64e{{\hslash }^{4}N}_0}{9{\left(2\pi \right)}^{1.5}x\left(1-x\right)E_{\mathrm{A}\mathrm{S}}^2{m_s}^{2.5}{\left(k_{\mathrm{B}}T\right)}^{0.5}}\psi \left(\eta \right)$$

(4)

$$\frac{1}{\mu }=\frac{1}{{\mu }_{\mathrm{A}\mathrm{P}}}+\frac{1}{{\mu }_{\mathrm{A}\mathrm{S}}}$$

(5)

$$\sigma =\mu en$$

(6)

$$\kappa ={\kappa }_{\mathrm{l}\mathrm{a}\mathrm{t}}+{\kappa }_{\mathrm{e}}={\kappa }_{\mathrm{l}\mathrm{a}\mathrm{t}}+L\sigma T$$

(7)

where $\mathrm{\hslash }$ is the reduced Planck constant,$F_i\left(\eta \right)$ the Fermi integral of order i, and $x$ is the alloy atomic composition in Sn, such as Mg₂Si_1−xSn_x ($x$ = 0.7 in this work), $\psi \left(\eta \right)=\frac{3\sqrt{\pi }}{16}\frac{F_{-0.5}\left(\eta \right)}{F_0\left(\eta \right)}$ and $L={\left(\frac{k_{\mathrm{B}}}{e}\right)}^2\frac{{3F}_0\left(\eta \right)F_2\left(\eta \right)-4{F_1}^2}{F_0{\left(\eta \right)}^2}$.

The other parameters are described in Table 1. These parameters were obtained from the literature for samples whose properties match ours [12], and are therefore applicable.

Figure 2. Comparison between the SPB modeled data (for n = 2.3∙10²⁶ m⁻³) based on the microscopic parameters in Table 1 and measured data: (a) Seebeck coefficient, (b) electrical conductivity, and (c) thermal conductivity.

Figure 2. Comparison between the SPB modeled data (for n = 2.3∙10²⁶ m⁻³) based on the microscopic parameters in Table 1 and measured data: (a) Seebeck coefficient, (b) electrical conductivity, and (c) thermal conductivity.

Table 1. Parameters used for the Single Parabolic Band model calculations of Mg₂Si_0.3Sn_0.7.

Parameter (Units)	Symbol	Value (SI)	Reference
Average density of states effective mass	$m_D^{*}$	2.1 $m_0$	[50]
Band degeneracy	$N_{\mathrm{v}}$	6	[31]
Single band mass	$m_s$	$m_D^{*}/N_v^{2/3}$	-
Theoretical mass density (g/cm3)	${\rho }_{\mathrm{D}}$	3.117	[50]
Longitudinal speed of sound (m/s)	$v_{\mathrm{l}}$	5290	linear with $x$, [50]
AP deformation potential constant (eV)	$E_{\mathrm{D}\mathrm{e}\mathrm{f}}$	9.8	[55]
Alloy scattering potential (eV)	$E_{\mathrm{A}\mathrm{S}}$	0.5	[31,55]
Number of atoms per unit volume (m−3)	$N_0$	4.105∙1028	linear with $x$, [31]

Table 1. Parameters used for the Single Parabolic Band model calculations of Mg₂Si_0.3Sn_0.7.

Parameter (Units)	Symbol	Value (SI)	Reference
Average density of states effective mass	$m_D^{*}$	2.1 $m_0$	[50]
Band degeneracy	$N_{\mathrm{v}}$	6	[31]
Single band mass	$m_s$	$m_D^{*}/N_v^{2/3}$	-
Theoretical mass density (g/cm3)	${\rho }_{\mathrm{D}}$	3.117	[50]
Longitudinal speed of sound (m/s)	$v_{\mathrm{l}}$	5290	linear with $x$, [50]
AP deformation potential constant (eV)	$E_{\mathrm{D}\mathrm{e}\mathrm{f}}$	9.8	[55]
Alloy scattering potential (eV)	$E_{\mathrm{A}\mathrm{S}}$	0.5	[31,55]
Number of atoms per unit volume (m−3)	$N_0$	4.105∙1028	linear with $x$, [31]

${\kappa }_{\mathrm{l}\mathrm{a}\mathrm{t}}$ is obtained from measured experimental data of our as-sintered n-type material: ${\kappa }_{\mathrm{l}\mathrm{a}\mathrm{t}}={\kappa }_{\mathrm{e}\mathrm{x}\mathrm{p}}-L{\sigma }_{\mathrm{e}\mathrm{x}\mathrm{p}}T\approx {\kappa }_{\mathrm{e}\mathrm{x}\mathrm{p}}-\left[1.5+\exp \left(-\frac{\left|{\alpha }_{\mathrm{e}\mathrm{x}\mathrm{p}}\right|}{116}\right)\right]{\sigma }_{\mathrm{e}\mathrm{x}\mathrm{p}}T$, where ${\kappa }_{\mathrm{e}\mathrm{x}\mathrm{p}}$, ${\sigma }_{\mathrm{e}\mathrm{x}\mathrm{p}}$, ${\alpha }_{\mathrm{e}\mathrm{x}\mathrm{p}}$ denote the bulk measured values, corresponding to the homogeneous material, before contacting and device making [11]. The obtained data is fitted with a third-order polynomial ${\kappa }_{\mathrm{l}\mathrm{a}\mathrm{t}}=1.63\left[\frac{\mathrm{W}}{\mathrm{m}\mathrm{K}}\right]-2.21·{10}^{-3}\left[\frac{\mathrm{W}}{\mathrm{m}{\mathrm{K}}^2}\right]·(T-273)+1.21·{10}^{-6}\left[\frac{\mathrm{W}}{\mathrm{m}{\mathrm{K}}^3}\right]·{(T-273)}^2+3.09·{10}^{-9}\left[\frac{\mathrm{W}}{\mathrm{m}{\mathrm{K}}^4}\right]·{(T-273)}^3$ with $T$ in K, and this fit equation is used to calculate the temperature- and position-dependent total thermal conductivity using the SPB model.

The comparison between bulk measurements on homogeneous samples corresponding to the n-type material employed for the TEG and the related SPB model results (for n = 2.3 ∙ 10²⁶ m⁻³, obtained with the experimental room temperature Seebeck coefficient and Equation (2)) and measured data, is shown in Figure 2.

It can be seen that SPB represents the experimental data well and captures the temperature dependence of the transport properties. $\alpha \left(T\right)$ starts to show a different temperature dependence only for $T\ge 400°\mathrm{C}$ due to the excitation of minority carriers; which, is beyond the analyzed temperature range. Therefore, SPB is fully valid for the initial, homogeneous material.

Obtaining the spatial carrier concentration profile:

The spatial variation of the Seebeck coefficient is obtained using an in-house-built device called the Potential & Seebeck Microprobe (PSM) [56,57]. This device locally measures the Seebeck coefficient and the voltage along a conductive sample; which, in the case of a TE leg, allows us to calculate the electrical contact resistance and map the Seebeck coefficient. Exemplary line scans, obtained in the PSM for the legs prior to TEG making, are shown in Figure 3. In [29], we measured an n-type leg of the TEG post cycling, and little difference was observed in the PSM Seebeck coefficient profile, showing a relatively stable behavior of the legs through TEG measurement. Therefore, pre-cycling profiles can be used for the SPB calculations in this work.

The Seebeck coefficient values obtained with the PSM are systematically underestimated, and not as accurate as those obtained under integral measurement conditions using the HTSσ device. This is due to the temperature difference between the effective position of the thermocouple junction and the point where the thermovoltage is measured [56]; which, leads to an empirically determined deviation between 10% and 20% of the measured Seebeck values in the PSM. This deviation depends on sample properties and tip wear, and was found for the range of thermal conductivities of TE materials (2 to 6.5 W/(mK), respectively, for Bi₂Te₃ and FeSi₂); in which, the range our inhomogeneous material falls. This effect is also known as the cold finger effect [43,51]. Note though, that the effect of statistical noise can be reduced by averaging.

A spatial Seebeck coefficient profile at room temperature is obtained by scanning an n-type leg with the PSM. The obtained profile ${\alpha }_{\mathrm{P}\mathrm{S}\mathrm{M}}\left(x\right)$ is converted into a “true” $\alpha \left(x\right)$ profile using previously measured Seebeck coefficient values from integral measurement conditions ${\alpha }_{\mathrm{i}\mathrm{n}\mathrm{t}}$. The employed assumptions here are first, that there is a constant relative difference between ${\alpha }_{\mathrm{P}\mathrm{S}\mathrm{M}}$ and $\alpha$ (which is plausible as the cold finger effect on the PSM measurements leads to a constant relative error) and second, that the carrier concentration in the middle of the sample (x = x₀) is unchanged due to the distance to outside metallization layers. This allows us to obtain a corrected Seebeck coefficient profile from $\alpha \left(x\right)={\alpha }_{\mathrm{P}\mathrm{S}\mathrm{M}}\left(x\right)\ast {\alpha }_{\mathrm{i}\mathrm{n}\mathrm{t}}/{\alpha }_{\mathrm{P}\mathrm{S}\mathrm{M}}\left(x_0\right)$, where ${\alpha }_{\mathrm{P}\mathrm{S}\mathrm{M}}\left(x_0\right)$ = −90 µV/K and ${\alpha }_{\mathrm{i}\mathrm{n}\mathrm{t}}$ = −109 µV/K. The deviation between those values is 17%; which, lies within the combined uncertainty for the local Seebeck coefficient by the PSM (10–20%) and that of the integral Seebeck coefficient measurement system (5%).

The assumption that the carrier concentration in the middle of the sample is unchanged is, in principle, an uncertainty of the CPIM. If that is not true, the correction factor and subsequent analysis are flawed. However, for the considered example, this assumption is reasonable since the PSM value at x₀ is close to values measured for as-sintered pellets made with other powder batches of the same composition and similar properties.

As described in Equations (1) and (2), the Seebeck coefficient at room temperature can be linked to the carrier concentration using the SPB model. A carrier concentration profile n(x) can therefore be obtained from the corrected Seebeck coefficient profile.

The temperature function T(x) is obtained along the leg assuming a linear profile between T_h,TE and T_c,TE (see CPM section). For Mg₂X materials, the linearity of the temperature profile can be assumed in spite of the interplay between Thomson heat, κ(T), and Joule heat, as shown by Ponnusamy et al. [42,58].

Constant Property Model (CPM):

The basics of the CPM are given in [39,40,41] while those applied to the TEG are derived in detail in [29]; however, the most relevant equations are given below and the relevant parameters are reported in Table 2:

$$∆T_{\mathrm{T}\mathrm{E},0}=\frac{U_{0,\mathrm{m}}}{N({\alpha }_{\mathrm{p}}-{\alpha }_{\mathrm{n}})}$$

(8)

$$T_{\mathrm{h},\mathrm{T}\mathrm{E},I}=T_{\mathrm{h},\mathrm{m}}-0.5\times {∆T}_{\mathrm{p}\mathrm{a}\mathrm{r},I}$$

(9)

$$T_{\mathrm{c},\mathrm{T}\mathrm{E},I}=T_{c,\mathrm{m}}+0.5\times {∆T}_{\mathrm{p}\mathrm{a}\mathrm{r},I}$$

(10)

$$\frac{Q_{\mathrm{o}\mathrm{p}\mathrm{t},\mathrm{m}}}{Q_{0,\mathrm{m}}}=\frac{{∆T}_{\mathrm{p}\mathrm{a}\mathrm{r},\mathrm{o}\mathrm{p}\mathrm{t}}}{{∆T}_{\mathrm{p}\mathrm{a}\mathrm{r},0}}=\frac{{∆T}_{\mathrm{m},\mathrm{o}\mathrm{p}\mathrm{t}}-{∆T}_{\mathrm{T}\mathrm{E},\mathrm{o}\mathrm{p}\mathrm{t}}}{{∆T}_{\mathrm{p}\mathrm{a}\mathrm{r},0}}$$

(11)

$$Q_I=K_{\mathrm{T}\mathrm{E}}{∆T}_{\mathrm{T}\mathrm{E},I}+I·N·\left({\alpha }_{\mathrm{p}}-{\alpha }_{\mathrm{n}}\right)T_{\mathrm{h},\mathrm{T}\mathrm{E},I}-\frac{1}{2}{I}^{2}R$$

(12)

$$I_{\mathrm{o}\mathrm{p}\mathrm{t}}=\frac{N{(\alpha }_{\mathrm{p}}-{\alpha }_{\mathrm{n}}){∆T}_{\mathrm{T}\mathrm{E},\mathrm{o}\mathrm{p}\mathrm{t}}}{2R}$$

(13)

$$R=R_{\mathrm{T}\mathrm{E}}+R_{\mathrm{c}}+R_{\mathrm{C}\mathrm{u}}=N\left[\frac{{\rho }_{\mathrm{p}}L}{A_{\mathrm{p}}}+\frac{{\rho }_{\mathrm{n}}L}{A_{\mathrm{n}}}+2r_{\mathrm{c}}\left(\frac{1}{A_{\mathrm{p}}}+\frac{1}{A_{\mathrm{n}}}\right)\right]+R_{\mathrm{C}\mathrm{u}}$$

(14)

$$K_{\mathrm{T}\mathrm{E}}=N\left(\frac{{\kappa }_{\mathrm{p}}{A}_{\mathrm{p}}}{L}+\frac{{\kappa }_{\mathrm{n}}{A}_{\mathrm{n}}}{L}\right)$$

(15)

$$P_{\mathrm{m}\mathrm{a}\mathrm{x}}=N\left(P_{\mathrm{n}}+P_{\mathrm{p}}\right)=\frac{{\left(N\left({\alpha }_{\mathrm{p}}-{\alpha }_{\mathrm{n}}\right){∆T}_{\mathrm{T}\mathrm{E},\mathrm{o}\mathrm{p}\mathrm{t}}\right)}^2}{4R}$$

(16)

$${\eta }_{\mathrm{m}\mathrm{a}\mathrm{x}}=\frac{T_{\mathrm{h},\mathrm{T}\mathrm{E},\mathrm{o}\mathrm{p}\mathrm{t}}-T_{\mathrm{c},\mathrm{T}\mathrm{E},\mathrm{o}\mathrm{p}\mathrm{t}}}{T_{\mathrm{h},\mathrm{T}\mathrm{E},\mathrm{o}\mathrm{p}\mathrm{t}}}\frac{\sqrt{1+Z{T}_{\mathrm{m}}}-1}{\frac{T_{\mathrm{c},\mathrm{T}\mathrm{E},\mathrm{o}\mathrm{p}\mathrm{t}}}{T_{\mathrm{h},\mathrm{T}\mathrm{E},\mathrm{o}\mathrm{p}\mathrm{t}}}+\sqrt{1+Z{T}_{\mathrm{m}}}}$$

(17)

$$Z=\frac{{N^2(\alpha }_{\mathrm{p}}-{\alpha }_{\mathrm{n}})2}{K_{\mathrm{T}\mathrm{E}}R}$$

(18)

where:

Table 2. Parameters of the CPM.

Symbol	Description
m (subscript)	Indicates measured value
$U_0$	Seebeck voltage
$T_{\mathrm{h},\mathrm{m},I}$, $T_{\mathrm{c},\mathrm{m},I}$	Temperature at the hot, cold block in TEG measurement at current I
$T_{\mathrm{h},\mathrm{T}\mathrm{E},I}$, $T_{\mathrm{c},\mathrm{T}\mathrm{E},I}$	Temperature at the hot, cold side of the TE legs at current I
${∆T}_{\mathrm{p}\mathrm{a}\mathrm{r},I}={∆T}_{\mathrm{m},I}-{∆T}_{\mathrm{T}\mathrm{E},I}$	Parasitic temperature loss at current I
$Q_I$	Heat flow at the hot side ($Q$in) at current I
$K_{\mathrm{T}\mathrm{E}}$	Thermal conductance of the TE legs
$R_{\mathrm{T}\mathrm{E}}$	Electrical resistance of the TE legs
$R_{\mathrm{c}}$	Electrical contact resistance
$r_{\mathrm{c}}$	Electrical contact resistivity
$N$	Number of leg pairs
L, $A$	Length, Cross-sectional area of TE element
$R_{\mathrm{C}\mathrm{u}}={\displaystyle \sum _i}\frac{L_i}{{\sigma }_{\mathrm{C}\mathrm{u}}{\left(T\right)A}_i}$	Resistance of the Cu bridges: sum of the resistances of all i pieces (varying geometries and temperatures). $L_i$,$A_i$ are the length, cross-sectional area of each Cu piece.
$I_{\mathrm{o}\mathrm{p}\mathrm{t},\mathrm{p}}$	Current at maximum power
$P_{\mathrm{m}\mathrm{a}\mathrm{x}}$	Maximum power output
${\eta }_{\mathrm{m}\mathrm{a}\mathrm{x}}$	Maximum efficiency

Table 2. Parameters of the CPM.

Symbol	Description
m (subscript)	Indicates measured value
$U_0$	Seebeck voltage
$T_{\mathrm{h},\mathrm{m},I}$, $T_{\mathrm{c},\mathrm{m},I}$	Temperature at the hot, cold block in TEG measurement at current I
$T_{\mathrm{h},\mathrm{T}\mathrm{E},I}$, $T_{\mathrm{c},\mathrm{T}\mathrm{E},I}$	Temperature at the hot, cold side of the TE legs at current I
${∆T}_{\mathrm{p}\mathrm{a}\mathrm{r},I}={∆T}_{\mathrm{m},I}-{∆T}_{\mathrm{T}\mathrm{E},I}$	Parasitic temperature loss at current I
$Q_I$	Heat flow at the hot side ($Q$in) at current I
$K_{\mathrm{T}\mathrm{E}}$	Thermal conductance of the TE legs
$R_{\mathrm{T}\mathrm{E}}$	Electrical resistance of the TE legs
$R_{\mathrm{c}}$	Electrical contact resistance
$r_{\mathrm{c}}$	Electrical contact resistivity
$N$	Number of leg pairs
L, $A$	Length, Cross-sectional area of TE element
$R_{\mathrm{C}\mathrm{u}}={\displaystyle \sum _i}\frac{L_i}{{\sigma }_{\mathrm{C}\mathrm{u}}{\left(T\right)A}_i}$	Resistance of the Cu bridges: sum of the resistances of all i pieces (varying geometries and temperatures). $L_i$,$A_i$ are the length, cross-sectional area of each Cu piece.
$I_{\mathrm{o}\mathrm{p}\mathrm{t},\mathrm{p}}$	Current at maximum power
$P_{\mathrm{m}\mathrm{a}\mathrm{x}}$	Maximum power output
${\eta }_{\mathrm{m}\mathrm{a}\mathrm{x}}$	Maximum efficiency

Note that all material properties ($\rho ,S,\kappa$) in Equations (8)–(18) are actually temperature averages, e.g., ${\alpha }_{\mathrm{p}}=\frac{1}{∆T}{\int }_{T_{\mathrm{c}}}^{T_{\mathrm{h}}}{\alpha }_{\mathrm{p}}\left(T\right)\mathrm{d}T$, ${\rho }_{\mathrm{p}}=\frac{1}{∆T}{\int }_{T_{\mathrm{c}}}^{T_{\mathrm{h}}}{\rho }_{\mathrm{p}}\left(T\right)\mathrm{d}T$ and ${\kappa }_{\mathrm{p}}={\left(\frac{1}{∆T}{\int }_{T_{\mathrm{c}}}^{T_{\mathrm{h}}}{{\kappa }_{\mathrm{p}}}^{-1}\left(T\right)\mathrm{d}T\right)}^{-1}$.

When modeling module properties, the resistance of the bridges is often neglected, [22,59] as in our previous work [29]. Here, $R_{\mathrm{C}\mathrm{u}}$ was considered explicitly in Equation (14), and represents 4% of the calculated total resistance. This not-so-small value arises despite a relatively large bridge thickness of 250 µm because of its relatively long length $L_i$ and the significantly reduced conductivity of the Cu at the hot side temperature [60]. Note also that we have taken the total bridge length; which, systematically overestimates the effective length. The dimensions of the hot side bridges are $L_{\mathrm{h}}$ = 12 mm and $A_{\mathrm{h}}$ = 6 mm; those of the shorter bridge on the cold side are $L_{\mathrm{c},\mathrm{s}}$ = 12 mm and $A_{\mathrm{c},\mathrm{s}}$= 5 mm; and those of the two longer bridges are $L_{\mathrm{c},\mathrm{l}}$ = 16 mm and $A_{\mathrm{c},\mathrm{l}}$ = 5 mm. $r_{\mathrm{c}}=4$ µΩcm² is used for all legs, as measured and reported in [29]. Similarly to [29], temperatures at $I_{\mathrm{o}\mathrm{p}\mathrm{t},\mathrm{P}}$ are used also for ${\eta }_{\mathrm{m}\mathrm{a}\mathrm{x}}$ calculations because the difference between $I_{\mathrm{o}\mathrm{p}\mathrm{t},\mathrm{P}}$ and $I_{\mathrm{o}\mathrm{p}\mathrm{t},\mathrm{\eta }}$ < 7%, so temperature conditions at both currents are similar.

3. Results

The Seebeck coefficient, electrical resistivity and thermal conductivity profiles of the n-type legs are obtained from $n\left(x\right)$ and $T\left(x\right)$. The procedure is represented in Figure 4.

From the calculated profiles, average material properties are calculated for the n-type leg; while for the p-type, a homogeneous material is assumed. The averages are obtained in dependence of the relevant temperature interval ($\mathrm{\Delta }{T}_{\mathrm{T}\mathrm{E},I}$), obtained using the CPM with the equations presented above, following, e.g., the procedure outlined in [29]. Module parameters are also calculated following those equations. A comparison of the relevant quantities assuming a homogeneous n-type leg and an inhomogeneous leg using the CPIM is presented in Figure 5.

Figure 4. Schematics of the procedure for the calculation of the property profiles for an inhomogeneous leg: (a) exemplary line scan of the Seebeck coefficient at room temperature, measured with PSM and also corrected for the cold finger effect; (b) carrier concentration (spatial) profile obtained from (a) using the SPB model; (c) fitted lattice thermal conductivity from experimental data of the sample directly after sintering; (d) corresponding Seebeck coefficient profile (in a temperature gradient) calculated using Equation (1); (e) corresponding electrical resistivity profile calculated using Equation (6); (f) corresponding thermal conductivity profile calculated using Equation (7). A comparative profile for a homogeneous material using the properties directly after sintering is added in dashed lines in (d–f). In (c–f), a linear temperature profile is assumed between T_c = 25 °C and T_h = 400 °C (for illustration purposes) to convert the temperature dependence of the properties into a spatial profile.

Figure 4. Schematics of the procedure for the calculation of the property profiles for an inhomogeneous leg: (a) exemplary line scan of the Seebeck coefficient at room temperature, measured with PSM and also corrected for the cold finger effect; (b) carrier concentration (spatial) profile obtained from (a) using the SPB model; (c) fitted lattice thermal conductivity from experimental data of the sample directly after sintering; (d) corresponding Seebeck coefficient profile (in a temperature gradient) calculated using Equation (1); (e) corresponding electrical resistivity profile calculated using Equation (6); (f) corresponding thermal conductivity profile calculated using Equation (7). A comparative profile for a homogeneous material using the properties directly after sintering is added in dashed lines in (d–f). In (c–f), a linear temperature profile is assumed between T_c = 25 °C and T_h = 400 °C (for illustration purposes) to convert the temperature dependence of the properties into a spatial profile.

Figure 5. Comparison of experimental data for a Mg₂(Si,Sn) TEG (measured during cooling sequence under stabilized temperature conditions, taken from [29]) to CPM and CPIM calculated data for a TEG with homogeneous and inhomogeneous n-type legs, respectively: (a) ratio of the temperature difference at the TE legs and the measured temperature difference at the module/heat flow meter interface at I = 0 and I_opt,P; (b) heat flows at I = 0 and I_opt,P; (c) inner resistance; (d) maximum conversion efficiency for measured inner resistance $R_m$ and several calculated inner resistances (considering TE resistance, contact resistance such as $r_{\mathrm{c}}=4$ µΩcm² and $R_{\mathrm{C}\mathrm{u}}$ ); (e) maximum power for measured and calculated inner resistances. The legend in (e) also applies in (d). $T_{\mathrm{h},\mathrm{m}}$ is the temperature at the hot side of the TEG in the measurement column.

Figure 5. Comparison of experimental data for a Mg₂(Si,Sn) TEG (measured during cooling sequence under stabilized temperature conditions, taken from [29]) to CPM and CPIM calculated data for a TEG with homogeneous and inhomogeneous n-type legs, respectively: (a) ratio of the temperature difference at the TE legs and the measured temperature difference at the module/heat flow meter interface at I = 0 and I_opt,P; (b) heat flows at I = 0 and I_opt,P; (c) inner resistance; (d) maximum conversion efficiency for measured inner resistance $R_m$ and several calculated inner resistances (considering TE resistance, contact resistance such as $r_{\mathrm{c}}=4$ µΩcm² and $R_{\mathrm{C}\mathrm{u}}$ ); (e) maximum power for measured and calculated inner resistances. The legend in (e) also applies in (d). $T_{\mathrm{h},\mathrm{m}}$ is the temperature at the hot side of the TEG in the measurement column.

Going from the homogeneous material to the inhomogeneous material, carrier loss near the metallization interfaces leads to an increase of the average absolute Seebeck coefficient and electrical resistivity while the thermal conductivity decreases. This trend is observed in the CPM results shown in Figure 5, where the temperature difference at the TE legs (based on measured open-loop voltage of the TEG and the input Seebeck coefficient) and the heat flows decrease for the inhomogeneous legs. Note that the significant decrease in heat flow (~10%) is due to a combined effect of the changed temperatures and the decreased $K$.

Employing Equation (14), the expected resistance of the TEG can be calculated from material properties, measured contact resistances before device making, and the contribution of the Cu bridges. If we consider the inhomogeneous leg properties the portion of the measured resistance that is due to the TE materials is 74%; while, it would be 70% if we consider the homogeneous properties of the original material. The contributions from contacts and bridge are visible, but relatively small. While the unexplained differences between measured and calculated resistances are large in both cases, this shows that incorrectly assuming homogeneous material would lead to an overestimation of these parasitic resistances. This corresponds to an overestimation of the “effective” contact resistances obtained from $r_{\mathrm{c},\mathrm{e}\mathrm{f}\mathrm{f}}=\frac{\left(R_{\mathrm{m}}-R_{\mathrm{l}\mathrm{e}\mathrm{g}\mathrm{s}}-R_{\mathrm{C}\mathrm{u}}\right)A_{\mathrm{n}}}{2N(1+A_{\mathrm{n}}/A_{\mathrm{p}})}$ [29] in the CPM model (at T_h = 200 °C: 65 µΩcm² for the homogeneous material, 58 µΩcm² for the inhomogeneous material). The relative difference in maximum power between inhomogeneous and homogeneous material (Figure 5e) is smaller for $R=R_{\mathrm{m}}$ than for $R=R_{\mathrm{i}\mathrm{d}\mathrm{e}\mathrm{a}\mathrm{l}}=R_{\mathrm{T}\mathrm{E}}+R_{\mathrm{c}}$ with $r_c=4$ µΩcm²; since for the latter, the resistance value is based on the calculated resistance of the TE materials, which differs; while for the first, the same measured resistance value is used for both materials.

The proportion of parasitic electrical resistance is large in both cases; this confirms that crack formation and propagation are responsible for the suboptimal performance of the TEG, as initially reported in [29]. By avoiding cracking, the maximum power output could be increased by 26%, with a corresponding power density of 1.06 W/cm². The smaller relative temperature difference for the inhomogeneous material also indicates that the thermal coupling between the heat source and TE material is worse than initially expected. From an application perspective, thermal losses inside the module and at the heat exchanger need to be minimized [61]. These are not reported for most module characterizations but will be relevant for system optimization and, hence, need to be characterized accurately.

It can be noted that the deviation from measured values for the maximum power in Figure 5e is similar for the homogeneous and for the inhomogeneous legs (2–4% and 3–5%, respectively, with increasing temperature); while, the inhomogeneous leg should be the most realistic case. This is due to the adjustment of the term $\left({\alpha }_{\mathrm{p}}-{\alpha }_{\mathrm{n}}\right){∆T}_{\mathrm{T}\mathrm{E},\mathrm{o}\mathrm{p}\mathrm{t}}$ to the measured $U_{0,\mathrm{m}}$ value; which, is common to both cases. It is expected that the CPIM error increases for higher temperatures where the SPB model is not as reliable, as the minority carriers start contributing to the conduction [62]. This limitation is even more pronounced for lower carrier concentrations where the conduction regime transition happens at a lower temperature. The hot-side portion of the leg is therefore where the SPB limitations and uncertainty are most relevant due to amplified bipolar effects; since, this portion combines higher temperature and lower carrier concentration. This will be further discussed below.

The calculated heat flow remains below the measured one and the deviation increases from 2–15% to 9–16% when going from a model considering a homogeneous to an inhomogeneous leg, due to the decrease of thermal conductivity with decreasing average carrier concentration of the material. In the case of inhomogeneous material, the deviation between the CPM-predicted heat flow and the measured value can therefore be occasionally larger than the estimated measurement uncertainty (13.5%). As explained in [63], the measurement uncertainty was obtained by testing commercial TEGs which have a larger number of legs, a higher filling factor, and a wider geometry [63]; therefore, it does not necessarily strictly apply to our TEG. As discussed in the SI of [29], the main challenge for small TEG prototypes is the radiative thermal bypass; which, could happen between the hot side and the cold side of the TEG itself, but also between the heater and the soldered cables and the heat flow meter (HFM) where the output heat flow is measured. Even though it cannot be stated with sufficient certainty, the improved analysis also indicates that the results of the efficiency measurements are systematically too low, i.e., the TEG performance is underestimated, as the otherwise quite predictive CPM model disagrees with the experimental data for heat flow and efficiency.

4. Discussion

Figure 6a shows the Seebeck coefficient with respect to temperature for experimental data and calculated data with the SPB model for different carrier concentrations. It can be seen that beyond a certain temperature the experimental data starts to bend much more than the SPB model due to the minority carrier contribution. This bending indicates the maximum temperature at which the SPB is reliable for each carrier concentration. The carrier concentration range shown in Figure 6a corresponds to the range determined in the inhomogeneous material, with lower carrier concentrations at the hot and cold sides and higher carrier concentrations in the middle portion. It can be therefore seen that, on an inhomogeneous leg with a gradient from 25 °C to 400 °C, the SPB prediction of the portion between approximately 350 °C and 400 °C would not be accurate.

Figure 6b,c show the calculated resulting profiles of the n-type material power factor (PF) and figure of merit ($zT$) assuming a linear temperature profile between 25 °C and 400 °C. As mentioned above, for a fraction of the leg on the hot side, the SPB data is likely not reliable; however, it can be seen that also the rest of the leg, the inhomogeneous leg is predicted to have a higher performance. This is only on a first glance in contradiction to the lower power CPIM calculation of the inhomogeneous leg in Figure 5e, as the reduced ${∆T}_{\mathrm{T}\mathrm{E}}$ for the inhomogeneous material overcompensates the increase in $\alpha$. Efficiency, on the other hand, is also governed by the heat flow, which is significantly lower in the inhomogeneous leg, leading to the increased efficiency with CPIM shown in Figure 5d.

Figure 6. (a) Seebeck coefficient with respect to temperature for: experimental data (from this work and the literature [64]) and calculated data for different carrier concentrations, obtained with the SPB model. It was added for visualization although no perfect match is expected as SPB parameters were obtained by different synthesis routes and for different sample compositions. The grey area visualizes the temperature range in which the SPB model apparently deviates from experimental data. Calculation of the (b) PF and (c) zT profiles for an inhomogeneous (SPB data) and a homogeneous leg (experimental data) with an assumed linear temperature profile between T_c = 25 °C and T_h = 400 °C.

Figure 6. (a) Seebeck coefficient with respect to temperature for: experimental data (from this work and the literature [64]) and calculated data for different carrier concentrations, obtained with the SPB model. It was added for visualization although no perfect match is expected as SPB parameters were obtained by different synthesis routes and for different sample compositions. The grey area visualizes the temperature range in which the SPB model apparently deviates from experimental data. Calculation of the (b) PF and (c) zT profiles for an inhomogeneous (SPB data) and a homogeneous leg (experimental data) with an assumed linear temperature profile between T_c = 25 °C and T_h = 400 °C.

The SPB predictions in Figure 6b,c show that the Zn-induced gradients initially present on the n-type legs after contacting are not detrimental to the TEG’s performance. The as-sintered material was not synthesized intentionally with a lower carrier concentration from the start because the SPB model tends to overestimate the figure of merit for low n (see, e.g., [50,65]). In practice, a lower-doped material might therefore have a lower performance. Also, the material was synthesized aiming for a maximized $z{T}_{\mathrm{m}\mathrm{a}\mathrm{x}}$; while $z{T}_{\mathrm{a}\mathrm{v}\mathrm{g}}$ is actually a better performance indicator [58].

In this work, a constant mobility parameter is assumed in the SPB model, independently of the carrier concentration value. If the local change in carrier concentration is due to Mg loss in the TE material (diffusion into the metallization), it could also have an effect on the carrier mobility, as was previously observed in [31,55]. For a one-dimensional current flow along the leg holds $\frac{\partial U}{\partial x}\propto \frac{1}{\sigma }$; an approximately constant mobility can therefore be verified by comparing the relative change in the slope of the electrical potential along a leg on a gradient portion and a middle portion of a leg as shown in Figure 3a in the light and dark grey zones (the slopes are, respectively, 8.3 mV/m and 3.5 mV/m) to the relative change of carrier concentration in those portions (8.5∙10²⁵ m⁻³ and 1.7∙10²⁶ m⁻³, respectively). In our case, both relative changes are about a factor of ≈2. This means that the electrical resistivity and the carrier concentration changed quite proportionally; therefore, mobility stayed relatively constant. This analysis validates the use of Equation (5) with a constant $E_{\mathrm{D}\mathrm{e}\mathrm{f}}$ for the inhomogeneous material and is of significant importance, as a change in mobility is one of the largest uncertainties of the SPB model.

Besides the limitations of the SPB model, part of the calculation error could also originate from the specific methodology of the CPIM, especially the determination of the corrected Seebeck coefficient profile, see the method section. Finally, the CPM is based on temperature averaging, which balances out the uncertainty on the carrier concentration profile; it is therefore, by definition, not sensitive to fine-tuning of the TE properties.

Lastly, the material change could have an impact on the self-compatibility of the material. This can be verified by calculating the compatibility factor $s\left(T\right)=\frac{(-1+\sqrt{1+zT})}{\alpha T}$ [38] for the minimum and maximum carrier concentration values of the profile; which, gives 4.4 for both carrier concentrations at 200 °C and 3.8 and 3.9 at 400 °C, respectively. Limited compatibility becomes an issue when $\mathit{s}$ differs by a factor of 2 [38]; it, therefore, is not expected to play a role here.

As TEGs are supposed to have a long-term application, a method of analysis that considers degradation mechanisms is crucial to understand device behavior. The new method presented here, CPIM, is easy to implement, as SPB-based material models are available for several relevant thermoelectric material systems [54,65,66] and the often-implemented CPM methodology needs to be adjusted only marginally. A further improvement of the presented approach would be to employ a two-band model [67] to overcome the limitations of the SPB model and Finite Element Modelling to cross-check and deepen the understanding of the results presented in this work.

5. Conclusions

In this paper, we present a new coupled material–device modeling approach, the CPIM, which is a methodology developed to capture property inhomogeneity in TE legs and which adapts the CPM model accordingly. This concept can be transferred to cases where the inhomogeneity is by design, for example for graded materials [38]. The CPIM furthermore enables an analysis of TEGs that have seen a thermal load where assuming homogeneous material could lead to wrong conclusions. For the considered case of a measured Mg₂(Si,Sn) based TEG, we obtain a smaller temperature difference at the TE legs and a reduced heat flow, compared to the initially assumed homogeneous material, beyond the assumptions and doubts of the SPB model. Consequently, our deepened analysis of the experimental data allows us to identify the heat flow measurement as a main future challenge for accurate measurement of small TEG prototypes.