The theory of thermoelectric effects has been studied extensively by many groups. With the advancement in material engineering and device fabrication, these devices have become more efficient and are being used in many applications. These devices are mainly divided into two major groups: The first group is based on Peltier effect     and is used in cooling and refrigeration applications   . These devices can be used in cooling or heating integrated circuits, biological and medical elements, medical sensors and devices, and other applications; The second group is based on Seebeck effect   and is mainly used in energy harvesting   -  . The energy harvesting devices can harvest wasted energy from factories, car exhausts, solar cells, power stations, and many other sources.
There are many types of materials that can be used to build thermoelectric devices and modules. These materials possess different ranges of Seebeck and Peltier effects  . Thermoelectric materials can be members of different groups. They can be from the Chalcogenides    , the Skutterudites  , the Half Heusler  , the Clathrates   , and the TAGS   . The most promising materials are the compounds based on Bismuth Telluride (Bi2Te3) and Antimony Telluride (Sb3Te3)   . Bi2Te3 and Sb3Te3 can be grown using many growth techniques such as: electrochemical  , electroplating  , flash evaporation  , thermal evaporation  and Pulsed Laser Deposition (PLD)     . Each technique has its advantages and disadvantages as discussed in  .
Thermoelectric devices have been modeled and simulated using several techniques. Spice program was used to simulate lumped and distributed parameter models   . Bipolar and diode thermoelectric devices have been modeled analytically for cooling and heat generation applications   . Thermoelectric devices have been modeled and simulated not only as individual devices but as a module for different applications   . These devices have been also simulated using numerical techniques  -  . Modeling and simulation are used to design devices, simulate their operation, and evaluate their performance before the costly fabrication process.
In this paper, we report on the analysis and simulation of multi-leg thermoelectric devices using Bi2Te3 and Sb3Te3. The study was conducted using the finite element (FE) multi-physics program, COMSOLTM  . The multi-leg devices under consideration were developed, fabricated, and tested by the authors  . In that study, electrical conductivity, Seebeck effect, and voltage sensitivity (SV) were measured  . The main goal of this study is to verify the accuracy of the experimental measurements and understand the device behavior with and without employing heat sinks. COMSOL provided the accuracy of a numerical analysis and the flexibility of modeling to emulate a real case scenario.
2. Numerical Model
The differential equations used in our analysis are part of COMSOL Program  . In this work, we used the thermoelectric module, where both heat and Poisson’s equations are used simultaneously to solve for the field variables temperature, T, and voltage, V,   . The partial differential equations used to obtain the device parameters are
where σ is the electrical conductivity, α is Seebeck coefficient, and k is the thermal conductivity of the material. The simulation model, built in COMSOL, is explained in extensive details in references   .
The structure of each device was constructed in COMSOL. The fabricated device structure and the properties of every material were used to build the model. Si (100) substrates, with a diameter of 3-inch and a thickness of 380 µm, were used. A3000 Å layer of SiO2was placed as insulator on top of the substrates. The thermoelectric materials, Bi2Te3 andSb2Te3, and the copper electrodes were added according to the device design. The thicknesses of these layers is taken to be400nm, which is their actual thickness in the fabricated devices  . Figure 1(a) shows a photograph of the fabricated devices. Figure 1(b) shows the device geometry constructed using COMSOL. Figure 2(a) shows the masks used in the fabrication, which indicate the areas of the Bi2Te3, Sb2Te3, and cupper layers. The top green-red-green and the bottom green-blue-green structures are for the characterization of Bi2Te3 and Sb2Te3, respectively. The cross structures are the registration marks for mask alignment. It should be noted that the test devices and the registration marks have not been added to the COMSOL geometry.
3. Simulation and Analysis of the Devices Tested without Heat Sink
Four devices, Device-1, Device-2, Device-3, and Device-4, whose thermoelectric materials were deposited at 100˚C, 200˚C, 300˚C and 400˚C, respectively, are simulated. The geometries were constructed with the actual device dimensions. The values of σ and k, used in the simulation analysis, are the measured values
Figure 1. (a) A photograph of the device after fabrication  ; and (b) The device geometry constructed using COMSOL.
Figure 2. The 8-leg thermoelectric device: (a) The masks for the three layers: the red is for the Bi2Te3, the blue is for Sb2Te3, and the green is for the copper metal; (b) Temperature distribution obtained from applying a 0.3 mV across terminals 1 and 2.
shown in Table 1. The values of Seebeck coefficient were assumed to be −200 ×10−6 V/K for Bi2Te3, 200 × 10−6 V/K for Sb2Te3, and 6.5 × 10−6 V/K for copper  . For copper electrodes, an electric conductivity of 5.9 × 108 and a thermal conductivity of 350 were used   . The thermal and electric conductivities for Bi2Te3 and Sb2Te3are obtained from measurements and published literatures, and they are listed in Table 1. All other physical materials were obtained from COMSOL library.
The simulation shows that when a voltage was applied across terminals I and 2, a temperature difference was generated between terminals 3 and 4, as shown in Figure 2(a). For device-1, when a voltage of 0.1 mV was applied, a temperature difference (ΔT) of 0.85 K was generated. And when the voltage increased to 0.2 mV and then to 0.3 mV, the temperature difference increased to 1.8 K and then to 2.67 K, respectively. Figure 2(b) shows one case of temperature simulation across the junctions, where an input voltage of 0.3 mV was applied. This voltage resulted in a temperature difference of 2.67 K. The figure shows also that the temperature difference is the same for each junction, while the electric potential decreases along the path from terminal 1 to terminal 2. This means that the junctions are thermally connected in parallel, but they are electrically connected in series. Simulations for other devices have produced similar results, but with different values. The values for the applied voltage and the corresponding temperature differences, for all devices, are shown in Table 2.
It has been observed from simulation that any increase in the applied voltage beyond certain values yielded no increase in the temperature difference. This appears to be a result of uncontrolled heat dissipation from the junctions, due to the absence of heat sink. The uncontrolled heat dissipation, made the device less sensitive to the applied voltage at higher values. This low sensitivity at higher voltages induced large errors, if SV or ST are obtained by curve fitting the entire V-T curve. Since the data at low voltages are more accurate, obtaining the
Table 1. Material properties used in the simulation [Ref. 24, 26, and 27 cited in  ].
Table 2. Temperature difference between junctions 3 and 4, when voltages were applied at terminals 1 and 2 of the devices obtained from COMSOL simulation.
sensitivities values by obtaining the slope of the curve at lower temperatures would yield more accurate values.
Similar to Device-1, for Device-2, applying any voltage above 0.4 mV did not result in a significant change in the temperature. For Device-3 and Device-4, 0.6 mV was the voltage at which no significant difference in temperature was observed. The exact behavior of the four devices was observed experimentally too. The agreement between experimental measurements and simulation data for all devices proves the validity of this study, and it is a strong evidence of the accuracy of our findings and conclusions.
To characterize the device performance, we simulated the voltage sensitivity (SV) and compared it to the corresponding experimental values. The sensitivity values for the simulation and experimental results are shown in Table 3. As shown in the table, the two sets of values are close to each other with about −7% to 20% difference. This deviation can be attributed to the use of constant values for Seebeck coefficient for all simulated devices; while in reality, Seebeck coefficient changes with the material growth conditions. The change of the Seebeck coefficient with growth temperature was automatically factored in the experimental measurements.
Temperature sensitivities (ST = 1/SV) for the devices were also calculated and their values are shown in Table 3, along with experimental results. As expected, the simulation values are in agreement with the experimental values, with deviation of −7% to 20%. Again this deviation can be attributed to the use of a
Table 3. Voltage Sensitivity (SV = ΔV/ΔT) and Temperature Sensitivity (ST = ΔT/ΔV) for experimental measurements  and COMSOL simulation (this work).
constant Seebeck coefficient for the thermoelectric materials, regardless of the growth temperature of the material. The decision to maximize SV versus maximizing ST depends on the application the device will be used for. For example, for power harvesting applications, SV should be maximized. On the other hand, for cooling/heating applications, ST should be optimized. In these cases, the device parameters, the geometries, and the fabrication process should be optimized for the intended goal.
4. Analysis of the Devices with Heat Sink
To obtain a better understanding of the performance, the devices were simulated again, but using a heat sink. In this case, the inner junctions, indicated as (I) in Figure 3, were kept at 293 K, while the temperature of the outer junctions, indicated as (O), varied from 293 K to 323 K. The electromotive force (emf), or the potential, created across terminals 1 and 2 was obtained. In addition, the heat flux across the junctions in both the x and y directions were simulated. The heat flux was simulated in order to understand the direction of heat flow, and to know whether heat is confined to the junctions or dissipate to the surroundings.
As can be seen from Table 4, the temperature of the outer junctions increased from T = 293 K to T = 323 K in a 6-degree increment. This resulted in temperature difference from 0 to 30 K. The temperature difference, in turn, generated a voltage difference across terminals 1 and 2, as shown in Figure 2. Table 4 shows the voltage across the two terminals, the temperature of the hot junctions, and the temperature difference between the hot and the cold junctions. To obtain the voltage sensitivities, the voltage generated across terminals 1 and 2 are plotted against the temperature, as shown in Figure 4.
As can be seen from Figure 4, the relationship between the applied voltage and the temperature is slightly nonlinear, bending to the right. However, the relationship can still be approximated as a straight line. The slopes of these lines represent the voltage sensitivities (SV) for the devices. Table 5 shows the device voltage sensitivities for Device-1, Device-2, Device-3, and Device-4 to be 175.56, 163.33, 177.77, and 187.77 µV/K, respectively. The Table shows also that the temperature sensitivities are 5.7, 6.12, 5.62, and 5.33 K/mV, respectively. Since ST is the reciprocal of SV, it is expected that devices with high voltage sensitivities will have low temperature sensitivities. The value of the sensitivity that matters
Figure 3. Thermoelectric device structure indicating the position of the inner and the outer parts of junctions.
Figure 4. Voltage generated across terminals 1 and 2 when temperature difference is maintained across the device junctions: (a) Device-1; (b) Device-2; (c) Device-3; and (d) Devices-4, with the temperature of the inner junctions is fixed at 293 K.
Table 4. Voltage generated across terminals 1 and 2 as a result of temperatures difference across the junctions. The temperature of the inner junctions is kept at 293 K.
depends on the application the device is intended for. For devices designed for cooling applications, the values of ST should be maximized. But for devices intended for energy harvesting, the value of SV should be the one to be maximized.
COMSOL simulation of Device-3 with outer junction temperature of 317 K is shown in Figure 5, as an example of the simulations results for this case. Similar results were obtained for all four devices, as explained before. Figure 5(a) shows the Temperature distribution, Figure 5(b) shows the Voltage distribution, Figure 5(c) shows the Heat Flux in the x-direction, and Figure 5(d) shows the Heat
Table 5. COMSOL simulation of Voltage Sensitivity (SV = ΔV/ΔT) and Temperature Sensitivity (ST = ΔT/ΔV) for devices with heat sink.
Figure 5. COMSOL simulation of Device-3 with outer junction temperature of 317 K. (a) Temperature distribution; (b) Voltage distribution (c); Heat Flux in the x-direction, and (d) Heat Flux in the y-direction.
Flux in the y-direction. The figure shows that the cold junction temperature is set at 293 K and the hot junction temperature is set at 317 K. From Figure 5(a), it can be seen that there is a temperature gradient across the junctions. This means that there is heat flowing from the hot to the cold junctions, generating electrical potential. The value of the generated potential is shown in Figure 5(b).
As can be seen from Figure 5(b), the voltage potential is the highest at terminal 1 and decreases until it reaches the lowest value at terminal 2. This decrease indicates that the junctions are electrically connected in series. By examining the potential values, we notice that the potential along the copper interconnects is constant. This means that there is little or no voltage drop across the metal. This means also that the emf, generated across terminals 1 and 2, is the summation of the voltage generated by all junctions, with no voltage drop across the series resistance.
To investigate heat dissipation or heat conductivity, the heat flux in the x-direction and the y-direction was simulated. From Figure 5(c), it is clear that there is heat flux across the junctions that are aliened in the x-direction only. Similarly, for the y-component, the flux is only in the junctions that are aliened along the y-direction; Figure 5(d). The Figure shows also that the heat flux in the substrate, between the individual legs, is zero. This indicates that heat exchange takes-place across the junction, and heat sinking is maintained.
All four devices that are analyzed here have been fabricated in our laboratory, but were tested without a heat sink. It was noticed that after certain voltage, voltage increase did not yield the expected temperature increase. It is believed that the lack of a heat sink allowed both junctions to change their temperature in the direction of thermal equilibrium. The exact phenomenon was also observed in the simulation results. Since the experimental and simulation results are in agreement, it can be claimed that both studies are accurate.
The voltage sensitivity, SV, and temperature sensitivity are used to characterize the device performance. For experimental results, SV increased from 93 µV/K to 132 µV/K for substrate temperature of 100˚C to 400˚C. On the other hand, for simulation, SV increased from 112 µV/K to 135 µV/K for substrate temperature of 100˚C to 400˚C. Although the values of SV for both cases are close, the increase of SV with growth temperature is sharper in case of experimental measurements than simulation. This mismatch can be attributed to the fact that simulation was performed assuming that Bi2T3 and Sb2T3 have constant Seebeck coefficients, regardless of the growth temperature of the materials. In reality, Seebeck coefficient changes with substrate temperature, which is reflected in the experimental measurements, but is not in the simulation analysis.
The simulation results show that sensitivity still increased with substrate temperature, even with using constant Seebeck coefficient. The increase could be due to the increase in the electrical conductivity and the decrease in the thermal conductivity of Bi2Te3 with growth temperature   . For Sb2Te3, the thermal conductivity is assumed constant and the electrical conductivity was found experimentally to decrease with the growth temperature, but with slower rate  . Even with these opposing effects, a moderate increase in the voltage sensitivity with growth temperature was observed.
When a heat sink is utilized, SV increases with growth temperature, as shown by simulation. SV increased from 175.56 to 187.77 µV/K. The increase of SV for devices with heat sink can be attributed to the absence of a load across terminals 1 and 2. Connecting a load results in drawing current, which, in turn, results in a voltage drop across the device series resistance. This, in turn, reduces the induced voltage across terminals 1 and 2, and consequently, reduces SV to values closer to those obtained from the experimental measurements. But since, in our simulation, no load was connected at the output, the induced voltage increased, resulting in a higher SV, as our results indicate.
By studying the sensitivities, one may realize that a device with low voltage sensitivity has high temperature sensitivity and assumes that such a device is most suitable for cooling and heating applications. This may not be true in reality, since device performance, for any application, depends on many factors, including, but not limited to, sensitivity. Hence, device structure, geometry, material properties, and processing should be designed to yield the best results for the intended application.
In conclusion, this study shows that reliable testing of thermoelectric devices requires the use of heat sink, to avoid uncontrolled change in the junction temperatures. The study shows also that the output load characteristics affect the device performance. Therefore, for a more accurate analysis of the devices, simulation and experimental measurement should be done with a load connected to the output terminals. The study has also shown that the device should be optimized for the application intended for. For energy harvesting, SV should be optimized. However, for cooling or heating, ST should be the parameter to be optimized. These parameters can be optimized by choosing the appropriate material and the optimum design structure that reduces stray elements.
Based on this study, two future works can be proposed: The first is to investigate the influence of device parameters on performance and the stability of the device operation. This will lead to developing design methodology that can be used to achieve the maximum intended performance; The second future work is to investigate the transient response and operation in the time domain. This would be essential for device design, if the device is used as a fast response heat sensor. Small signal analysis requires the development of a small signal equivalent circuit.
The authors would like to acknowledge Northern Illinois University and its Microelectronic Research and Development Laboratory (MRDL) for supporting this work.
 Abdel-Motaleb, I. and Qadri, S.M. (2017) Fabrication and Characterization of PLD-Grown Bismuth Telluride (Bi2Te3) and Antimony Telluride (Sb2Te3). Thermoelectric Devices. Journal of Electronics Cooling and Thermal Control, 7, 63-77.
 Tritt, T.M. (2002) Thermoelectric Materials: Principles, Structure, Properties, and Applications. Encyclopedia of Materials: Science and Technology, 2nd Edition, 1-11.
 Trinidad, P.M.P. and Carbajal, G. (2015) Potential Use of Thermoelectric Generator Device for Air Conditioning System. Proceedings of the 13th Latin American and Caribbean Conference for Engineering and Technology, Santo Domingo, 29-31 July 2015.
 Weerasinghe, R. and Hughes, T. (2017) Numerical and Experimental Investigation of Thermoelectric Cooling in Down-Hole Measuring Tools: A Case Study. Case Studies in Thermal Engineering, 10, 44-53.
 LeBlanc, S., Yee, S.K., Scullin, M.L., Dames, C. and Goodson, K.E. (2014) Material and Manufacturing Cost 564 Considerations for Thermoelectric. Renewable and Sustainable Energy Reviews, 32, 313-327.
 Chen, L.G., Meng, F.K. and Sun, F.R. (2016) Thermodynamic Analyses and Optimization for Thermoelectric Devices: The State of the Arts. Science China-Technological Sciences, 59, 442-455.
 Seetawan, T., Singsoog, K., Srichai, S, Thanachayanont, C., Amornkitbamrung, V. and Chindaprasirt, P. Thermoelectric Energy Conversion of p-Ca3Co4O9/n-CaMnO3 Module. Energy Procedia, 61, 1067-1070.
 Ma, Q., Fang, H., and Zhang, M. (2017) Theoretical Analysis and Design Optimization of Thermoelectric Generator. Applied Thermal Engineering, 127, 758-764.
 Cao, Q., Luan, W. and Wang, T. (2017) Performance Enhancement of Heat Pipes Assisted Thermoelectric Generator for Automobile Exhaust Heat Recovery. Applied Thermal Engineering, in Press.
 Chung, D.-Y., Hogan T., Brazis P., Rocci-Lane, M., Kannewurf, C., Bastea, M., Uher, C. and Kanatzidis, M.-G. (2000) CsBi4Te6: A High-Performance Thermoelectric Material for Low-Temperature Applications. Science, 287, 1024-1027.
 Uher, C., Yang, J., Hu, S., Morelli, D.T. and Meisner, G.P. (1999) Transport Properties of Pure and Doped MNiSn, M = (Zr, Hf). Physical Review B, 59, 8615-8621.
 Nolas, G.S., Cohn, J.L., Slack, G.A. and Schujman, S.B. (1998) Semiconducting Ge Clathrates: Promising Candidates for Thermoelectric Applications. Applied Physics Letters, 73, 178.
 Meng, J.F., Chandra Shekar, N.V., Badding, J.V. and Nolas, G.S. (2001) Threefold Enhancement of the Thermoelectric Figure of Merit for Pressure Tuned Sr8Ga16Ge30. Journal of Applied Physics, 89. 1730.
 Levin, E.M., Bud’ko, S.L. and Schmidt-Rohr, K. (2012) Enhancement of Thermopower of TAGS-85 High-Performance Thermoelectric Material by Doping with the Rare Earth Dy. Advanced Functional Materials, 22, 2766-2774.
 Lim, S.K., Kim, M.Y. and Oh, T.S. (2009) Thermoelectric Properties of the Bismuth-Antimony-Telluride and the Antimony-Telluride Films Processed by Electrode Position for Micro-Device Applications. Thin Solid Films, 517, 4199-4203.
 Kimi, M.-Y. and Oh, T.-S. (2011) Processing and Thermoelectric Performance Characterization of Thin-Film Devices Consisting of Electrodeposited Bismuth Telluride and Antimony Telluride Thin-Film Legs. Journal of Electronic Materials, 40, 759-764.
 Foucaran, A., Sackda, A., Giani, A., Pascal-Delannoy, F. and Boyer, A. (1998) Flash Evaporated Layers of (Bi2Te3-Bi2Se3)(N) and (Bi2Te3-Bi2Se3)(P). Materials Science and Engineering: B, 52, 154-161.
 Carmo, J.P., Goncalves, L.M., Wolffenbuttel, R.F. and Correia, J.H. (2010) A Planar Thermoelectric Power Generator for Integration in Wearable Microsystems. Sensors and Actuators A: Physical, 161, 199-204.
 Bailini, A., Donati, F., Zamboni, M., Russo, V., Passoni, M., Casari, C.S., Li Bassi, A. and Bottani, C.E. (2007) Pulsed Laser Deposition of Bi2Te3 Thermoelectric Films. Applied Surface Science, 254, 1249-1254.
 Mitrani, D., Salazar, J., Turó, A., García, M.J. and Chávez, J.A. (2007) Lumped and Distributed Parameter SPICE Models of TE Devices Considering Temperature Dependent Material Properties. 13th International Workshop on Thermal Investigation of ICs and Systems (THERMINIC), Budapest, Hungary, September 2007, 17-19.
 Chen, M., Rosendahl, L.A., Condra, T.J. and Pedersen, J.K. (2009) Numerical Modeling of Thermoelectric Generators with Varying Material Properties in a Circuit Simulator. IEEE Transactions on Energy Conversion, 24, 112-124.
 Ortega, P. and Olivares-Robles, M. (2017) Analysis of a Hybrid Thermoelectric Microcooler: Thomson Heat and Geometric Optimization. Entropy, 19, 312.
 Headings, L.M. (2011) Modeling and Development of Thermoelectric Device Technologies for Novel Mechanical Systems. PhD Dissertation, Ohio State University, Columbus, OH.
 Karami Lakeh, H.K., Abardeh, R.H. and Kaatuzian, H. (2015) Numerical Simulation of a Segmented Thermoelectric Generator. Energy and Power, 5, 1-9.
 Yan, D., Dawson, F., Pugh, M. and El-Deib, A. (2014) Time Dependent Finite Volume Model of Thermoelectric Devices. IEEE Transaction on Industry Applications, 50, 600-608.
 Seetawan, T., Seetawan, U., Ratchasina, A., Srichaia, S., Singsoog, K., Namhongsa, W., Ruttanapun, C. and Siridejachai, S. (2012) Analysis of Thermoelectric Generator by Finite Element Method. Procedia Engineering, 32, 1006-1011.
 COMSOL Multi-Physics.
 Zou, H., Rowe, D.M. and Williams, S.G.K. (2002) Peltier Effect in a Co-Evaporated (Sb2Te3)(P)-Bi2Te3(N) Thin Film Thermocouple. Thin Solid Films, 408, 270-274.