An Extensive Unified Thermo-Electric Module Characterization Method

Thermo-Electric Modules (TEMs) are being increasingly used in power generation as a valid alternative to batteries, providing autonomy to sensor nodes or entire Wireless Sensor Networks, especially for energy harvesting applications. Often, manufacturers provide some essential parameters under determined conditions, like for example, maximum temperature difference between the surfaces of the TEM or for maximum heat absorption, but in many cases, a TEM-based system is operated under the best conditions only for a fraction of the time, thus, when dynamic working conditions occur, the performance estimation of TEMs is crucial to determine their actual efficiency. The focus of this work is on using a novel procedure to estimate the parameters of both the electrical and thermal equivalent model and investigate their relationship with the operating temperature and the temperature gradient. The novelty of the method consists in the use of a simple test configuration to stimulate the modules and simultaneously acquire electrical and thermal data to obtain all parameters in a single test. Two different current profiles are proposed as possible stimuli, which use depends on the available test instrumentation, and relative performance are compared both quantitatively and qualitatively, in terms of standard deviation and estimation uncertainty. Obtained results, besides agreeing with both technical literature and a further estimation method based on module specifications, also provides the designer a detailed description of the module behavior, useful to simulate its performance in different scenarios.


Introduction
TEMs are increasingly being used for both power generation and for cooling/heating applications. In the first case the Seebeck effect is predominant to convert temperature gradient ∆T into electrical energy, whereas, in the second case, the dual effect-the Peltier effect-is exploited. TEMs' several positive features make them reliable [1,2], and thus their use is appealing for both academic and industrial researchers.
In particular, Thermo-Electric Generators (TEGs) can generate clean energy in a small space and therefore are being used in an increased number of standalone applications ranging from aerospace [3,4], to industry [5][6][7]. Their possible integration with other energy harvesting technologies was also investigated; for example, for increasing the overall efficiency of photovoltaics [8] by exploiting the thermal gradient between the back of the panels and the ambient temperature, hence recovering the otherwise dissipated heat. The inherent restrained necessity of maintenance and the high long-term efficiency, make TEGs a valid alternative to batteries in providing an autonomous source of energy to sensor nodes [9][10][11] or entire Wireless Sensor Networks [12][13][14][15][16].
In addition, Thermo-Electric Coolers (TECs) have been considered, for a long time, for a number of applications such in that described in [17], where they are used to keep the chip temperature below an assigned threshold in electronic packaging systems, or in specialized air conditioning [18,19] and refrigeration applications [7,[20][21][22].
Performances of a TEC are typically considered in correspondence to an optimal operating current that yields the maximum heat flux for a given operating condition, but in many cases, a TEC-based system is operated under the best conditions for only a fraction of the time, unless it is operated at a fixed heat load and fixed ambient conditions. Thus, the performance of a TEC at off-peak conditions are of great importance in the design phase and in determining the actual efficiency under dynamic working conditions. However, manufacturers often only provide some essential parameters to characterize TEMs. Regarding the specifications of a TEC, some parameters are specified in correspondence of a couple of temperature values of the hot side (typically 27 • C and 50 • C). These usually are the largest temperature differential ∆T max that can be obtained between the two sides along with the input direct current I max and the voltage U max that have produced it. Conversely, for TEG modules they usually provide the power P e m , the voltage U m and the maximum efficiency η m at matched load, i.e., when the load resistance R m equals the internal electrical resistance of the module. Clearly, these values do not correspond to actual ones but to maximum theoretical specifications, scarcely significant for practical use in many applications.
The need to quantitatively characterize TEMs' behavior, implies the design of a robust testbed capable of addressing a twofold objective: to estimate the parameters of a TEM for design purposes and to define a robust method for comparing different TEM models. The development of a useful testbed for the characterization of TEMs is the main topic of this work. The focus is on using an appropriate procedure and suitable hardware to investigate the underlying relationship between both the operating temperature, the temperature difference and the parameters of the electrical and thermal equivalent model for commercial modules. The novelty consists in applying a proper current stimulus (for which two alternatives are discussed) to the TEM under test and simultaneously acquire both thermal and electrical data to automatically derive all parameters in a single test, without using either a discrete load [9,23,24] nor hot plates or auxiliary TEMs [2,11,17,[24][25][26][27][28][29]. Unlike other approaches, the developed testbed permit, using a very simple configuration, to quickly obtain all relevant parameters, making possible an extensive TEM characterization over a wide range of temperature differences, ambient temperatures and electric loads. This paper is structured as follows: after the introduction of the thermoelectric Figure of Merit Z, a typical method to estimate the thermal resistance which involves the measurement of four different temperatures between the surfaces of a reference material is described in Section 2.1. The state-of-the art review about electrical characterization is then described in Section 2.2 and subsequently the proposed method is briefly explained in Section 2.3. In Sections 3 and 4 the electrical and thermal characterization are analytically approached, respectively, whereas in Section 5 the testbed is presented and detailed. Finally, results obtained with the two alternatives of the proposed approach are experimentally compared to each other in Section 6, confirming was expected from the analytical analysis. The obtained parameters are presented and validated with a third estimate based on the manufacturers' specification [30]. Section 7 summarizes the paper's content drawing the conclusions.

State of the Art and Proposed Method
The thermoelectric figure of merit Z is a parameter which summarizes the bulk material proprieties and allows the comparison among different TEMs. Alternatively, the same meaning is given to a dimensionless parameter ZT where T is the average absolute temperature between the two surfaces of the module. It is defined as: Sensors 2016, 16, 2114 3 of 20 Expression (1) requires the knowledge of the relative Seebeck coefficient α s , the internal electrical resistance R I N and the thermal resistance Θ I N of a TEM, all parameters that provide the most comprehensive description of the physical behavior of a TEM both as TEG and TEC. In fact, by recourse to them it is for instance possible to express the electrical power P e delivered to the load and the thermal power P t absorbed by the module in variable working conditions for a TEG module: where V L is the voltage on the load. Furthermore, the efficiency η of a module, i.e., the amount of thermal power converted in electrical power, may be derived in all working conditions as: If Equations (2) and (4) are calculated at match load conditions (i.e., R L = R I N ), it is possible to obtain the maximum value of the electrical power deliverable to the load P e m (or equivalently the maximum power density if dividing it by the module area) and the maximum efficiency η m which are generally reported in datasheets.
Methods used nowadays to estimate such parameters suffer from various limitations, such as the necessity to realize different experiments to carry out all the requested measurements for Z calculation [31], or conversely, they directly provide it without letting one extract the single components [32]. Anyway, none of these is able to provide an extensive description of TEMs behavior in variable working conditions. In the following, a summary of these methods is reported.

Thermal Characterization
Among the above mentioned three parameters, the thermal resistance Θ I N results, in principle, quite simple to be obtained but less immediate to be measured, since it requires different temperature probes and the realization of a suitable conditioning circuitry. It can be indirectly measured using a test configuration like the one described in Section 5, measuring the emitted and absorbed heat fluxes (q em and q abs , respectively) and applying the thermal equivalent of Ohm's law Θ = ∆T/q.
Once the fluxes are known, it is possible to write down the energy balance (5) between the Peltier effect, the heat conduction and Joule effect due to the injected current: where Θ re f is the thermal resistance of the reference medium. Its value can be previously fine-tuned sourcing the module with a known electrical power and solving the energy balancing equation q em − q abs = V I with no thermal losses considered since they are negligible if the measurement setup is thermally isolated. Of course, it would be necessary to know the thermal dependence of Θ re f for the adopted material with the mean temperature T. However, if the investigated temperature range is not too wide, a material with an almost constant thermal conductivity in the given range may be adopted. Therefore, in order to measure Θ I N , it is necessary to measure the current sourced to the module and both the emitted and the absorbed heat flux, which in turn require four temperatures to be measured or equivalently two voltage difference induced across a reference medium with a known thermal resistance Θ re f . In Equation (5), T h and T c refer to those temperatures directly measured at the ends of the internal thermoelectric elements (pellets) and can be computed taking into account the temperature drops on the ceramic surfaces of a TEM: where Θ cer is the thermal resistance of the ceramic surfaces. However, in many cases, the use of alumina oxide (Al 2 O 3 ) ceramic external plates in the module fabrication, makes the correction (6) negligible, since the joint combination of a high thermal conductivity (one order of magnitude higher than Bi 2 Te 3 ) and a small thickness (one order of magnitude lower than Bi 2 Te 3 pellets) lead to a thermal resistance Θ cer value that is usually negligible as furtherly demonstrated in Equation (20). Finally, solving Equation (5) with respect to Θ I N , it is possible to derive a measurement formula for the unknown thermal resistance: where α s is supposed to have been previously estimated, whereas I st is the steady-state current injected into the module under test. This approach is very common and widely adopted when thermal characterization is required, using simple or more complex setups, depending on the application and uncertainty requirements. For example, in [24,25] the TEM is placed between a heat source given by a controlled hot-plate and a controlled cooling fan, whereas in [26] a shielded heating block and a cold-plate are used to assess the TEM efficiency, monitoring the heat flux, the delivered power and also the mechanical load to make a good thermal contact between each component. Another interesting setup is proposed in [27], where both the heat produced by an electrical resistance and the cooling capacity provided by a TEM, are controlled by a proportional-integral controller. Other propose the use of one or two opportunely driven auxiliary TEMs to set the temperature gradient on the TEM under test [2,28,33]. For example, in [28], Θ I N is obtained measuring only one heat flux induced into an aluminum block with a known thermal conductivity, interposed between the module under test and an auxiliary TEM used to set the average temperature, hence requiring only three temperatures to be measured.

Electrical Characterization
At the end of fifties, Harman [34] developed the method, hence the name "Harman Method" (HM), used to obtain the constituent parameters of the thermoelectric figure of merit using both AC and DC measurements; the former to obtain the electrical resistivity whereas the latter to measure the Seebeck coefficient and the thermal conductivity. Generally, the HM is valid under the hypothesis of thermal equilibrium: In particular, when the temperature of the hot side T h , the cold side T c and the ambient one T a do not deviate from the average temperature T by more than 0.3%, ensuring a good linear approximation of the radiative component in heat transmission, accurate to 1% [34]. This condition is typically achieved by injecting a current, large enough to be measured but sufficiently small (generally in the order of milliamperes) to meet the condition (8). In this condition it can be proved [34] that the figure of merit ZT can be written as ZT = V α /V r , where V r is the ohmic component of the voltage across the sample and V α is the Seebeck voltage. In the testbed proposed by Harman, these two voltage components were alternatively measured immediately after switching a mechanical chopper on and off, responsible for the creation of a square-wave current. The main useful feature introduced by HM, was indeed the capability to estimate directly ZT in terms of a ratio of two voltages however, because of the technology limitations at the time, it suffered from a lack of precision and reproducibility [35]. The method proposed by Harman was further developed by Buist in the so called "Transient Method" (TM) [35], in which the TEM is driven by a pulsed current like the one reported in Figure 1. At steady-state, the current is switched on, resulting in an immediate increase in the measured voltage because of the occurrence of the ohmic voltage component V r . However, due to the slower characteristic response of heat phenomena compared to electrical ones, a temperature difference arises across the sample, generating a Seebeck voltage that slowly increases V r to V r + V α . voltages however, because of the technology limitations at the time, it suffered from a lack of precision and reproducibility [35]. The method proposed by Harman was further developed by Buist in the so called "Transient Method" (TM) [35], in which the TEM is driven by a pulsed current like the one reported in Figure 1. At steady-state, the current is switched on, resulting in an immediate increase in the measured voltage because of the occurrence of the ohmic voltage component . However, due to the slower characteristic response of heat phenomena compared to electrical ones, a temperature difference arises across the sample, generating a Seebeck voltage that slowly increases to + . The main difference between the HM and the TM is in the use of a modern instrumentation such as a DAQ board and a PC to reconstruct the entire voltage waveform over time ( Figure 1). The use of the TM implies that, for each electrical working condition, a different current pulse have to be generated and acquired, but also that to have a good repeatability, a sufficient number of voltage waveform have to be processed. The success of the TM is widely demonstrated by the large scientific literature concerning the TEM characterization that was inspired to it. Only to cite some: in [32], a test facility called "Z-meter" is proposed to measure the performance parameters of a TEM and some correction factors were introduced to take into account real environment conditions; in [28], the TM is used to estimate and in a commercial TEM. Another improvement with respect to the HM was proposed in [29], where both the temperature difference induced by applying a DC current and the heat flow created by a low frequency AC current, are used to obtain , and , whereas Θ is indirectly estimated using Equation (1). In [30] the HM is used to obtain the parameters of the TEM equivalent PSPICE model and to validate a method for extracting such model parameters from specifications. This latter feature is also used by the authors to derive a rough estimate of the TEM parameters and furtherly comparing them with the values obtained using the method proposed in this work.

Proposed Method
To overcome the limitations of the methods cited above, we propose a procedure which relies on two different modified versions of the TM for the electrical part while the typical approach described in Section 2.1 is used for the thermal characterization. The proposed method provides a detailed description of the three parameters that compose the figure of merit with respect to variations of the operating conditions, i.e., the temperature of both sides of the module.
Briefly, the proposed method consists in placing the TEM module under test between two heat-flux sensors; when the desired operating conditions are met, a small signal or a rapid current variation is applied and the following quantities are measured: the emitted and the absorbed fluxes, the temperature at cold and hot side ( and ), the voltage and the current at its terminals. The developed testbed automatically sweeps along a wide range of drive currents using only a The main difference between the HM and the TM is in the use of a modern instrumentation such as a DAQ board and a PC to reconstruct the entire voltage waveform over time ( Figure 1). The use of the TM implies that, for each electrical working condition, a different current pulse have to be generated and acquired, but also that to have a good repeatability, a sufficient number of voltage waveform have to be processed.
The success of the TM is widely demonstrated by the large scientific literature concerning the TEM characterization that was inspired to it. Only to cite some: in [32], a test facility called "Z-meter" is proposed to measure the performance parameters of a TEM and some correction factors were introduced to take into account real environment conditions; in [28], the TM is used to estimate α s and R I N in a commercial TEM. Another improvement with respect to the HM was proposed in [29], where both the temperature difference induced by applying a DC current and the heat flow created by a low frequency AC current, are used to obtain α s , R I N and ZT, whereas Θ I N is indirectly estimated using Equation (1). In [30] the HM is used to obtain the parameters of the TEM equivalent PSPICE model and to validate a method for extracting such model parameters from specifications. This latter feature is also used by the authors to derive a rough estimate of the TEM parameters and furtherly comparing them with the values obtained using the method proposed in this work.

Proposed Method
To overcome the limitations of the methods cited above, we propose a procedure which relies on two different modified versions of the TM for the electrical part while the typical approach described in Section 2.1 is used for the thermal characterization. The proposed method provides a detailed description of the three parameters that compose the figure of merit ZT with respect to variations of the operating conditions, i.e., the temperature of both sides of the module.
Briefly, the proposed method consists in placing the TEM module under test between two heat-flux sensors; when the desired operating conditions are met, a small signal or a rapid current variation is applied and the following quantities are measured: the emitted and the absorbed fluxes, the temperature at cold and hot side (T h and T c ), the voltage V L and the current I at its terminals. The developed testbed automatically sweeps along a wide range of drive currents using only a data-acquisition (DAQ) board and a 4-quadrant transconductance amplifier, but if the latter is not available, similar results may be obtained using the small signal method. Although all above parameters are measured at the same time, in the following section they are discussed separately to avoid confusion between electrical and thermal ones.

Electrical Measurement
A TEM can be electrically modelled as a voltage source V th = α s ∆T in series with an internal resistance R I N . Once the TEM reaches a steady-state condition, the voltage V L at its terminals is measured along with the current I flowing in it using a shunt resistor R s placed in series (described in Section 5). In this condition, the total voltage resulting at the power supply output terminals is: If a rapid current variation is applied, fast enough to not interfere with its thermal steady-state, then V th can be considered constant in such time window and the total voltage (9) is dependent only by the supplied current, hence on the electrical response of the module. Therefore, the resulting voltage at the terminals of the module under test is: To separate the two effects, this condition is obtained when the current variation occurs in a time window τ el which is at least one order of magnitude smaller than the thermal time constant τ th .

Current Sweep (CS) Method
The adopted current variation consists in a CS from the value I st it has in the steady-state condition up to its opposite −I st . Using Equation (10), this CS produces a line drawn in the I − V plane, as shown in Figure 2, where the slope is given by R I N and V th is its intercept with the ordinate axis, both referred to a given temperature T c and T h . data-acquisition (DAQ) board and a 4-quadrant transconductance amplifier, but if the latter is not available, similar results may be obtained using the small signal method. Although all above parameters are measured at the same time, in the following section they are discussed separately to avoid confusion between electrical and thermal ones.

Electrical Measurement
A TEM can be electrically modelled as a voltage source = Δ in series with an internal resistance . Once the TEM reaches a steady-state condition, the voltage at its terminals is measured along with the current flowing in it using a shunt resistor placed in series (described in Section 5). In this condition, the total voltage resulting at the power supply output terminals is: If a rapid current variation is applied, fast enough to not interfere with its thermal steady-state, then can be considered constant in such time window and the total voltage (9) is dependent only by the supplied current, hence on the electrical response of the module. Therefore, the resulting voltage at the terminals of the module under test is: To separate the two effects, this condition is obtained when the current variation occurs in a time window which is at least one order of magnitude smaller than the thermal time constant .

Current Sweep (CS) Method
The adopted current variation consists in a CS from the value it has in the steady-state condition up to its opposite − . Using Equation (10), this CS produces a line drawn in the − plane, as shown in Figure 2, where the slope is given by and is its intercept with the ordinate axis, both referred to a given temperature and .

Small Signal (SS) Method
As described further in Section 4, the consistence of the CS method has been validated comparing it with a second current profile that produces similar results, but consists of a small sinusoidal signal added to the steady-state current , which can be considered to all effects as a perturbation around the bias point. From the theoretical point of view, such a small current does not require to satisfy condition (10) on the thermal time constant, since its amplitude is generally negligible with respect to .

Small Signal (SS) Method
As described further in Section 4, the consistence of the CS method has been validated comparing it with a second current profile that produces similar results, but consists of a small sinusoidal signal i ss added to the steady-state current I st , which can be considered to all effects as a perturbation around the bias point. From the theoretical point of view, such a small current does not require to satisfy condition (10) on the thermal time constant, since its amplitude is generally negligible with respect to I st .

Analytical Comparison among CS and SS Method
Intuitively, the SS method would seem a better approach since it does not require any additional hypothesis, but the choice of a CS as current profile provides far more uniform results and a higher reliability in the estimation of both the involved parameters. This statement can be easily demonstrated analytically. Indeed, Equation (10) is a linear regression problem [36] where a given set of n pairs The regression slope R IN , the intercept V th and respective standard deviations are computed using the following expressions [37]: where I denotes the mean value of I, SSI = ∑ n i=1 I i − I 2 is the sum of I squared errors, s 2 I is the sample, not population, variance of I and s 2 V L ,I is the error variance. From Equation (14), it is immediate to see that the standard deviation σ V th in the estimation of the intercept V th for a given σ R IN , increases linearly together with the samples standard deviation s I and with the mean value I, which is always zero if a symmetric CS is used. Conversely, I = I st when using a SS current profile and this produces higher errors when the value of the steady-state current is increased. Similarly, the variance σ 2 R IN in the estimation of the slope R IN is inversely dependent on the variance of I and increases as long as the samples tend to form a spot instead of spreading on the plane. This therefore means that using a symmetric CS, hence acquiring wider spread samples, leads to reduced errors in the estimation of fit parameters with respect to the SS method, but also that the CS method has a decreasing uncertainty in the slope estimation for increasing currents, whereas the SS method has a constant uncertainty.
After computing the standard deviations σ V th and σ R IN , if the estimation errors are unbiased, normally and independently distributed, the 100(1 − α)% confidence intervals on the slope and intercept can be expressed considering that the errors in the estimations of the same are both distributed as t random variables [38] with n − 2 degrees of freedom: If the applied CS is not sufficiently fast to satisfy Equation (10), both because of an erroneous over-estimation of the module thermal time constant (Section 4) or because of intrinsic limits in the acquisition speed of the implemented test, a correction may be applied to the acquired signals. In fact, as shown further in Section 6, a quadratic term compares in the Kirchhoff's voltage law equation, that induce a bending of the I − V line while the current varies. In this case, all the acquired values can still be used for parameters estimation, but the actual slope is the one at the beginning of the sweep, before the bending appears.
Hence, if Equation (10) is weakly satisfied, a linear regression with pure quadratic formulâ V Li = a + bI i + cI 2 i may be applied and the slope R IN and intercept V th may be computed as follows: where:

Thermal Characterization
Thermal measurements have been carried out using two heat flux sensors while varying the temperature of the two sides sourcing the module with an adequate current ( Figure 3). Each one was implemented by the authors, as described in a previous work [38], stacking two metal layers (2 mm thick) for temperature sensing with an interposed medium with known thermal conductivity k re f , both with the same geometric area of the module under test. In each metal layer, a small hole was made to accommodate the thermocouple sensors, whereas the contact resistance between layers was minimized using a thin layer of high thermal conductivity silver-based thermal paste. In order to satisfy Equation (10), the thermal transfer function of the module have to be preliminarily estimated; in this work it was used an Autoregressive Exogenous method included in the MATLAB ® System Identification Toolbox [39].
After applying a current step I(t), four temperatures (T 1 , T 2 , T 3 , T 4 ) are measured and composed to create two differential temperature signals: where ∆T 32 should be most representative of the response of the TEM itself, whereas ∆T 41 should also take into account the response of entire system of which it is part. Equations (15) are then used to estimate two transfer function, expressed as the Laplace transform of a first order system: where Aτ th is the DC gain.
The output differential temperature time constant at external sides of the heat-flux sensors is always greater than the one at direct contact with the TEM because of the heat propagating time in a given medium (i.e., τ th 41 τ th 32 ). Thus, τ th 32 is taken into consideration to establish the maximum temporal duration of the electrical perturbation τ el .
During the CS, all the four temperatures are acquired and elaborated according to the procedure described in Section 2.2 to compute the internal thermal resistance of the module Θ I N .
Since, even in the worst examined case (i.e., in correspondence of the highest heat fluxes), the temperature drops on the ceramic plates of the module (20) are negligible, it results that that T c ≈ T 2 and T h ≈ T 3 as previously stated in Equation (6): Making also explicit in Equation (4) the dependence of the remaining parameters on the acquired temperatures, we can rewrite the thermal resistance equation as a function of the measured variables as: where k re f = 1/Θ re f is the thermal conductivity of the reference medium, obtained as described previously in Section 2.1 as calibration factor, minimizing the squared error associated to the flux/power balancing equation. Therefore, it is necessary to adopt a material that exhibits a low dependence of its thermal conductivity with the temperature in the investigated range. At this point, if the following simplifying assumptions are made: • k re f has no associated uncertainty since it is a calibration factor • the four temperatures have the same uncertainty u T because they are measured with the same kind of temperature sensor (type-J thermocouple), they are acquired using the same DAQ board and referred to the same cold-junction sensor • all variables are uncorrelated, except I st and V th that are obviously correlated The uncertainty u Θ I N on the estimation of the internal thermal resistance can be derived applying standard uncertainty propagation [40] to Equation (21), hence, the thermal resistance uncertainty may be expressed as: from which is mainly observable a clear quadratic inverted trend with respect to ∆T 32 .

Automatic Test
The implemented automatic test procedure performs the above described measurement over a matrix of testing points obtained from the combination of two vectors containing respectively the desired temperatures for the cold side T c and the desired working ∆T. The whole test is acquired temperatures, we can rewrite the thermal resistance equation as a function of the measured variables as: where = 1 Θ ⁄ is the thermal conductivity of the reference medium, obtained as described previously in Section 2.1 as calibration factor, minimizing the squared error associated to the flux/power balancing equation. Therefore, it is necessary to adopt a material that exhibits a low dependence of its thermal conductivity with the temperature in the investigated range.
At this point, if the following simplifying assumptions are made:  has no associated uncertainty since it is a calibration factor  the four temperatures have the same uncertainty because they are measured with the same kind of temperature sensor (type-J thermocouple), they are acquired using the same DAQ board and referred to the same cold-junction sensor  all variables are uncorrelated, except and that are obviously correlated The uncertainty on the estimation of the internal thermal resistance can be derived applying standard uncertainty propagation [40] to Equation (21), hence, the thermal resistance uncertainty may be expressed as: from which is mainly observable a clear quadratic inverted trend with respect to  32 T .

Automatic Test
The implemented automatic test procedure performs the above described measurement over a matrix of testing points obtained from the combination of two vectors containing respectively the desired temperatures for the cold side and the desired working Δ . The whole test is conducted inside a Discovery Es 250 (DY-250) climate chamber (  All measured variables, outlined below, are voltage signals acquired at = 160 Hz by means of a 16-bit DAQ board X Series USB 6361 by National Instruments (Austin, TX, USA); in particular:


is the voltage at the TEM terminals  is the voltage across the shunt resistor All measured variables, outlined below, are voltage signals acquired at f s = 160 Hz by means of a 16-bit DAQ board X Series USB 6361 by National Instruments (Austin, TX, USA); in particular: • V L is the voltage at the TEM terminals • V s is the voltage across the shunt resistor R S • V a is the climate chamber ambient temperature measured by means of a LM35 temperature sensor (Texas Instruments, Dallas, TX, USA) • V 1 , V 2 , V 3 , V 4 are the voltages generated by the four different J-type thermocouples each one inserted into one surface of the heat-flux sensors • V cj is the cold junction compensation temperature acquired by means of a further LM35 sensor placed close to the DAQ board The same DAQ board is also responsible for generating the driving voltage for a 4-quadrant TOE7621 transconductance amplifier (TOELLNER Electronic Instrumente GmbH, Herdecke, Germany) that supplies the TEM under test.
The entire acquisition/generation process is handled by a LabVIEW Virtual Instrument (VI) that runs a proportional-integrative (PI) controller ( Figure 4) in a closed-loop feedback to drive the module to the desired steady state working conditions. Once known the current-to-thermal transfer function H ∆T (s) estimated in Equation (19), the tuning parameters for the PI controller have been computed using both the closed-loop Ziegler-Nichols method and the MATLAB ® Simulink Auto-tune tool. Although both methods produce good results, the latter has been chosen because of its ability to produce a desired response in a smaller time with less overshoot.


is the climate chamber ambient temperature measured by means of a LM35 temperature sensor (Texas Instruments, Dallas, TX, USA)  , , , are the voltages generated by the four different J-type thermocouples each one inserted into one surface of the heat-flux sensors  is the cold junction compensation temperature acquired by means of a further LM35 sensor placed close to the DAQ board The same DAQ board is also responsible for generating the driving voltage for a 4-quadrant TOE7621 transconductance amplifier (TOELLNER Electronic Instrumente GmbH, Herdecke, Germany) that supplies the TEM under test.
The entire acquisition/generation process is handled by a LabVIEW Virtual Instrument (VI) that runs a proportional-integrative (PI) controller ( Figure 4) in a closed-loop feedback to drive the module to the desired steady state working conditions. Once known the current-to-thermal transfer function estimated in Equation (19), the tuning parameters for the PI controller have been computed using both the closed-loop Ziegler-Nichols method and the MATLAB ® Simulink Auto-tune tool. Although both methods produce good results, the latter has been chosen because of its ability to produce a desired response in a smaller time with less overshoot. The steady-state conditions for a given combination of and Δ are met when the following two conditions are satisfied over 100 s or 16000 samples: where and are the mean and standard deviation of the measured temperature gradient, whereas and are the respective threshold values. In this case, the VI locks the driving current to the last current value, waits for a time equal to five thermal time constant ( = 5 ) to allow all transients to run out and then switches from Control Mode to Measure Mode ( Figure 5).
The measurement procedure consists of three sections run in sequence ( Figure 6), in which three different current profile are generated, each one with a duration that is set to a value sufficiently small with respect to the module thermal time constant ( ≤ ≪ ): (1) steady-state: the last driving current value is applied and the acquired values are used to compute statistics over . (2) SS: a small current 10 Hz sinusoidal stimulus is added to , with an amplitude equal to the ratio of standard deviation threshold to the static gain of Equation (19): (3) CS: the driving current is swept to its symmetric value − with a ramp-like signal that is finally switched back to the initial value. The steady-state conditions for a given combination of T c and ∆T are met when the following two conditions are satisfied over 100 s or 16,000 samples: where µ ∆T and σ ∆T are the mean and standard deviation of the measured temperature gradient, whereas µ th and σ th are the respective threshold values. In this case, the VI locks the driving current to the last current value, waits for a time equal to five thermal time constant (τ w = 5τ th 41 ) to allow all transients to run out and then switches from Control Mode to Measure Mode ( Figure 5). The measurement procedure consists of three sections run in sequence ( Figure 6), in which three different current profile are generated, each one with a duration τ m that is set to a value sufficiently small with respect to the module thermal time constant (τ m ≤ τ el τ th ): (1) steady-state: the last driving current value is applied and the acquired values are used to compute statistics over I st . (2) SS: a small current 10 Hz sinusoidal stimulus is added to I st , with an amplitude equal to the ratio of standard deviation threshold to the static gain of Equation (19): (3) CS: the driving current I st is swept to its symmetric value −I st with a ramp-like signal that is finally switched back to the initial value.  At the end of each test, all raw data (acquired voltage signals) are stored for further post-processing using MATLAB: is divided by to obtain .  and are obtained from respective voltages using the LM35 nominal sensitivity S = 10 mV/(°C).  , , , are computed using the coefficients provided by National Institute of Standards and Technology (NIST) and applying a software cold-junction compensation.  and are computed as described in Section 2.1. Once all the raw signals have been properly scaled, the three parameters composing the thermoelectric figure of merit can be calculated for each working condition , Δ :  At the end of each test, all raw data (acquired voltage signals) are stored for further post-processing using MATLAB: is divided by to obtain .  and are obtained from respective voltages using the LM35 nominal sensitivity S = 10 mV/(°C).  , , , are computed using the coefficients provided by National Institute of Standards and Technology (NIST) and applying a software cold-junction compensation.  and are computed as described in Section 2.1. Once all the raw signals have been properly scaled, the three parameters composing the thermoelectric figure of merit can be calculated for each working condition , Δ : At the end of each test, all raw data (acquired voltage signals) are stored for further post-processing using MATLAB: T a and T cj are obtained from respective voltages using the LM35 nominal sensitivity S = 10 mV/( • C).
• T 1 , T 2 , T 3 , T 4 are computed using the coefficients provided by National Institute of Standards and Technology (NIST) and applying a software cold-junction compensation. Once all the raw signals have been properly scaled, the three parameters composing the thermoelectric figure of merit ZT can be calculated for each working condition {T c , ∆T}: • V th and R IN are derived using Equation (10), applying a linear regression to acquired values {I, V L }. • α s is then obtained as ratio of V th to ∆T. • Θ I N is computed using Equation (7).

Experimental Results
The above described procedure has been used to characterize the performance of a commercial low-cost module TES1-12730 from Thermonamic Electronics Corporation (Nanchang, China) [41]. For sake of clarity, the experimental results are split in sub-sections: firstly, the results of the module identification procedure will be presented, then the consistence of the SS method is shown and finally the CS method is validated and characterization results will be shown.

Module Identification
The TEM has been identified applying a 0.1 V voltage step to the amplifier, which results in a 320 mA current step because of the transconductance gain G m = 3.2 A/V. The measured differential output signals ∆T 32 and ∆T 41 are thus used to obtain the two respective transfer functions, which continuous-time models are resumed in Table 1; H 32 ( Figure 7) is then used to tune the feedback loop in Figure 4.  and are derived using Equation (10), applying a linear regression to acquired values , .  is then obtained as ratio of to Δ .  Θ is computed using Equation (7).

Experimental Results
The above described procedure has been used to characterize the performance of a commercial low-cost module TES1-12730 from Thermonamic Electronics Corporation (Nanchang, China) [41]. For sake of clarity, the experimental results are split in sub-sections: firstly, the results of the module identification procedure will be presented, then the consistence of the SS method is shown and finally the CS method is validated and characterization results will be shown.

Module Identification
The TEM has been identified applying a 0.1 V voltage step to the amplifier, which results in a 320 mA current step because of the transconductance gain = 3.2 A/V. The measured differential output signals Δ and Δ are thus used to obtain the two respective transfer functions, which continuous-time models are resumed in Table 1; (Figure 7) is then used to tune the feedback loop in Figure 4.

SS Method Analysis
After elaborating the acquired data for the SS current stimulus, statistics have shown the consistence of this measurement method. The linear regression model associated with Equation (10), for each of the 210 tests, was estimated with a mean value of the determination coefficient Before starting the test, the adopted value of the time test window τ m = 1 s has been considered sufficiently small compared with the thermal time constant τ th 32 = 11.5 s of H 32 . This indeed have produced N s = τ m f s = 160 samples for each measurement section.

SS Method Analysis
After elaborating the acquired data for the SS current stimulus, statistics have shown the consistence of this measurement method. The linear regression model associated with Equation (10), for each of the 210 tests, was estimated with a mean value of the determination coefficient R 2 = 0.99550 and a standard deviation σ R 2 = 0.00023. If we consider the test with the worst regression fit (R 2 = 0.9919) relative to T C = 49 • C and ∆T = 45 • C, the plot of residuals (Figure 8) shows that they are, in any case, normally distributed with zero mean.  Furthermore, the probability plot in Figure 9a, shows how the distribution of the estimation residuals compared with a normal distribution with matched variance, confirms a reasonable fit to normally distributed residuals with almost no evident outliers. The residual lag plot (Figure 9b) also shows the absence of any remaining serial correlation among residuals, since they are all distributed uniformly among the four quadrants. Hence, the SS method seems to satisfy all demands for the application of the linear regression and the confidence intervals relative to the estimated parameters may be expressed as in Equation (15). In Figure 10, the and standard deviations obtained using Equations (13) and (14) are reported for all working conditions, showing that the is almost constant through all tests. In fact, as explained in Section 3.3, it mainly depends on the variances of data, which is set to a constant value in Equation (24). Similarly, the linear trend in the intercept standard deviation is due to its dependence on the mean value of data, which increases for increasing Δ . Furthermore, the probability plot in Figure 9a, shows how the distribution of the estimation residuals compared with a normal distribution with matched variance, confirms a reasonable fit to normally distributed residuals with almost no evident outliers. The residual lag plot (Figure 9b) also shows the absence of any remaining serial correlation among residuals, since they are all distributed uniformly among the four quadrants.  Furthermore, the probability plot in Figure 9a, shows how the distribution of the estimation residuals compared with a normal distribution with matched variance, confirms a reasonable fit to normally distributed residuals with almost no evident outliers. The residual lag plot (Figure 9b) also shows the absence of any remaining serial correlation among residuals, since they are all distributed uniformly among the four quadrants. Hence, the SS method seems to satisfy all demands for the application of the linear regression and the confidence intervals relative to the estimated parameters may be expressed as in Equation (15). In Figure 10, the and standard deviations obtained using Equations (13) and (14) are reported for all working conditions, showing that the is almost constant through all tests. In fact, as explained in Section 3.3, it mainly depends on the variances of data, which is set to a constant value in Equation (24). Similarly, the linear trend in the intercept standard deviation is due to its dependence on the mean value of data, which increases for increasing Δ . Hence, the SS method seems to satisfy all demands for the application of the linear regression and the confidence intervals relative to the estimated parameters may be expressed as in Equation (15). In Figure 10, the R I N and V th standard deviations obtained using Equations (13) and (14) are reported for all working conditions, showing that the σ R I N is almost constant through all tests. In fact, as explained in Section 3.3, it mainly depends on the variances s 2 I of I data, which is set to a constant value in Equation (24). Similarly, the linear trend in the intercept standard deviation is due to its dependence on the mean value of I data, which increases for increasing ∆T.

CS Method Analysis
Applying to CS data the same procedure just described, the obtained main determination coefficient is = 0.999965 with a standard deviation = 2.93 × 10 .
Despite of the high value of , the residuals of the linear regression for the same working point analyzed before (i.e., = 49 °C and Δ = 45 °C), show a non-zero mean and a clear non-normal distribution as evinced from Figure 11a, which means that a linear fit is not able to fully explain the underlying data relation. This is probably due to the fact that the adopted measure time window , even if one order of magnitude smaller then thermal time constant, is too large to satisfy the hypothesis in Equation (10). Nevertheless, a visual inspection of residuals plot (Figure 11b) shows a clear quadratic relation, which suggests to fit the data using a non-linear formula.
Hence, applying the least squares method with a quadratic formula as described at the end of Section 3.3, better fit results are obtained ( Figure 11a) with a resulting mean determination coefficient = 0.999996 and a standard deviation = 2.22 × 10 . Also, the resulting residuals mean error is zero and these are normally distributed, confirming a successful fit.
(a) (b) Figure 11. (a) Normal probability plot of the residuals of the fitted model: linear (blue) and quadratic (red); (b) Leverage plot of data and model As mentioned in Section 3.3, the CS method produces far more uniform results with much lower uncertainties, as showed in Figure 12.

CS Method Analysis
Applying to CS data the same procedure just described, the obtained main determination coefficient is R 2 = 0.999965 with a standard deviation σ R 2 = 2.93 × 10 −6 .
Despite of the high value of R 2 , the residuals of the linear regression for the same working point analyzed before (i.e., T c = 49 • C and ∆T = 45 • C), show a non-zero mean and a clear non-normal distribution as evinced from Figure 11a, which means that a linear fit is not able to fully explain the underlying data relation. This is probably due to the fact that the adopted measure time window τ m , even if one order of magnitude smaller then thermal time constant, is too large to satisfy the hypothesis in Equation (10). Nevertheless, a visual inspection of residuals plot ( Figure 11b) shows a clear quadratic relation, which suggests to fit the data using a non-linear formula.

CS Method Analysis
Applying to CS data the same procedure just described, the obtained main determination coefficient is = 0.999965 with a standard deviation = 2.93 × 10 .
Despite of the high value of , the residuals of the linear regression for the same working point analyzed before (i.e., = 49 °C and Δ = 45 °C), show a non-zero mean and a clear non-normal distribution as evinced from Figure 11a, which means that a linear fit is not able to fully explain the underlying data relation. This is probably due to the fact that the adopted measure time window , even if one order of magnitude smaller then thermal time constant, is too large to satisfy the hypothesis in Equation (10). Nevertheless, a visual inspection of residuals plot (Figure 11b) shows a clear quadratic relation, which suggests to fit the data using a non-linear formula.
Hence, applying the least squares method with a quadratic formula as described at the end of Section 3.3, better fit results are obtained ( Figure 11a) with a resulting mean determination coefficient = 0.999996 and a standard deviation = 2.22 × 10 . Also, the resulting residuals mean error is zero and these are normally distributed, confirming a successful fit.
(a) (b) Figure 11. (a) Normal probability plot of the residuals of the fitted model: linear (blue) and quadratic (red); (b) Leverage plot of data and model As mentioned in Section 3.3, the CS method produces far more uniform results with much lower uncertainties, as showed in Figure 12. Hence, applying the least squares method with a quadratic formula as described at the end of Section 3.3, better fit results are obtained (Figure 11a) with a resulting mean determination coefficient R 2 = 0.999996 and a standard deviation σ R 2 = 2.22 × 10 −6 . Also, the resulting residuals mean error is zero and these are normally distributed, confirming a successful fit. As mentioned in Section 3.3, the CS method produces far more uniform results with much lower uncertainties, as showed in Figure 12. Since when using this method, the variance of is not constant and increases for higher Δ , a clear downward trend occurs in ; whereas , which depends both on ̅ and , shows a growing trend because of the increment of the latter, while the former tends to zero due to the sweep symmetry.
Conversely, the standard deviation of the thermal resistance ( Figure 13) follows a trend similarly to as expected from Equation (22).

Experimental Comparison among CS and SS Method
In order to test the compatibility among the CS and SS method, for all working conditions, a point by point relative error Equation (25) is computed and mean values as well as standard deviations are reported in  (25) where xx P represent the generic estimated TEM parameter (i.e., , , Θ or ). Since when using this method, the variance of I is not constant and increases for higher ∆T, a clear downward trend occurs in σ R I N ; whereas σ V th , which depends both on I and s 2 I , shows a growing trend because of the increment of the latter, while the former tends to zero due to the sweep symmetry.
Conversely, the standard deviation of the thermal resistance σ Θ I N ( Figure 13) follows a trend similarly to σ R I N as expected from Equation (22). Since when using this method, the variance of is not constant and increases for higher Δ , a clear downward trend occurs in ; whereas , which depends both on ̅ and , shows a growing trend because of the increment of the latter, while the former tends to zero due to the sweep symmetry.
Conversely, the standard deviation of the thermal resistance ( Figure 13) follows a trend similarly to as expected from Equation (22).

Experimental Comparison among CS and SS Method
In order to test the compatibility among the CS and SS method, for all working conditions, a point by point relative error Equation (25) is computed and mean values as well as standard deviations are reported in  (25) where xx P represent the generic estimated TEM parameter (i.e., , , Θ or ).

Experimental Comparison among CS and SS Method
In order to test the compatibility among the CS and SS method, for all working conditions, a point by point relative error Equation (25) is computed and mean values as well as standard deviations are reported in Table 2: e r = P CS −P SS P SS e r = 1 n i n j ∑ i∈T c ∑ j∈∆T P CSij −P SSij P SSij (25) where P xx represent the generic estimated TEM parameter (i.e., α S , R I N , Θ I N or ZT). Computed mean relative errors show a good match between the two methods, with halved error values for CS method using a quadratic formula with respect to SS method and a linear regression.
The result of the characterization process of the module TES1-12730, using the proposed CS method, produced interesting results. Both the Seebeck coefficient α S (Figure 14a) and the internal electrical resistance R I N (Figure 14b) increase their value almost linearly as long as T C ∝ T and ∆T increase. Despite the Seebeck voltage appears smooth and uniform, the trend of the Seebeck coefficient shows to be more irregular because it is derived as ratio of the former to the temperature difference ∆T, which is less accurate for small values.  Computed mean relative errors show a good match between the two methods, with halved error values for CS method using a quadratic formula with respect to SS method and a linear regression.
The result of the characterization process of the module TES1-12730, using the proposed CS method, produced interesting results. Both the Seebeck coefficient ( Figure 14a) and the internal electrical resistance (Figure 14b) increase their value almost linearly as long as ∝ and Δ increase. Despite the Seebeck voltage appears smooth and uniform, the trend of the Seebeck coefficient shows to be more irregular because it is derived as ratio of the former to the temperature difference Δ , which is less accurate for small values. Regarding the internal thermal resistance Θ (Figure 15a), it shows a strong dependence with increasing Δ , whereas it seems to be only slightly affected by variations of , which decrease induces an observed weak resistance increment. The roughness of the surface for small Δ is due also in this case to its analytical dependence on measured temperatures. Regarding the internal thermal resistance Θ I N (Figure 15a), it shows a strong dependence with increasing ∆T, whereas it seems to be only slightly affected by variations of T C , which decrease induces an observed weak resistance increment. The roughness of the surface for small ∆T is due also in this case to its analytical dependence on measured temperatures. Regarding the internal thermal resistance Θ (Figure 15a), it shows a strong dependence with increasing Δ , whereas it seems to be only slightly affected by variations of , which decrease induces an observed weak resistance increment. The roughness of the surface for small Δ is due also in this case to its analytical dependence on measured temperatures. The figure of merit ZT (Figure 15b) exhibits a strong linear dependence on ∆T and a weaker one on T C , and since it is a combination of the previous mentioned parameters, it also shows a bending that is taken into account using second order terms in the interpolating function.

Conclusions
The novelty of the proposed technique consists in using a simple testbed to stimulate the module and to simultaneously acquire electrical and thermal data so to obtain the three parameters that compose the figure of merit Z. An extensive characterization of TEMs in many variable working conditions is hence performed within a single test.
The developed testbed automatically sweeps along a wide range of drive currents using only a DAQ board and a 4-quadrant transconductance amplifier, but if the latter instrument is not available, the authors have shown that similar results may also be obtained using a different approach.
To validate the method, we firstly compared the results obtained using two different current profiles: a small signal one that require no previous assumption and a sweep signal that conversely requires a preliminary hypothesis to be satisfied, but that produce far more accurate results; furtherly we compared them with those provided in the datasheet by the manufacturer. The performance of both approaches was theoretically analyzed and experimentally verified, demonstrating a good match between the two methods through a comparison among data generated by CS and SS profiles.
If using the current sweep approach, the authors have shown that obtained results may be furtherly improved without redoing all the tests, but simply applying a numerical correction to acquired data.
The experimental results (i.e., α S , R I N , Θ I N and ZT) for variable ∆T and T C working conditions have been reported in exhaustive tri-dimensional graphs, which values (reported in Table 3) are in agreement with those extrapolated from the datasheet using a different method.
Finally, simple analytical relations were derived from interpolation of obtained data and reported in Table 4 to provide a simple tool for eventual comparison with other data and methods and to simulate the TEM performance in different scenarios.
The reported modeling at different operating conditions, may play an important role in estimate the performance of TEMs for energy generation, affecting for instance the harvestable energy as long as the internal electrical resistance moves away from designed matched value.

Conflicts of Interest:
The authors declare no conflict of interest.