Spatial Equivalent Circuit Model for Simulation of On-Chip Thermoelectric Harvesters

Interest in autonomous low-power energy sources has risen with the development and widespread use of devices with very low energy consumption. Interest in thermoelectric harvesters has increased against this background. Thermoelectric harvesters, especially harvesters on-chip, have peculiar properties related to the thermal route, thermal transients, and spatial temperature distribution within the chip. A behavioral model of the harvester is required for engineers to successfully develop voltage converters with maximum power point tracking and energy storage units. There are accurate models based on the finite element method, but these models are usually not compatible with simulators of electrical circuits, and therefore are not convenient for designers. Existing equivalent circuit models fit this requirement, but usually do not consider many parameters. This article proposes an original method that allows simulating spatial thermoelectric processes by analogy with the finite difference method, using electrical circuits simulations software. The study proposes a complete methodology for building the model and examples of simulations of one-, two- and three-dimensional problems, as well as examples of simulation of macro problems in the presence of external thermal and electrical devices, such as heatsink and electrical load.


Introduction
The thermoelectric properties of materials have been known, and used with greater or lesser intensity, for over a hundred years. Thermoelectric power generators steadily occupy a niche in the market for autonomous power supplies used in space [1,2], in the ocean [3][4][5] and in other hard-to-reach places. Thermoelectric harvesters that turn waste heat into energy have been further developed these days. Such harvesters are increasingly used as an autonomous power source in the nodes of wireless sensor networks. The newest trends in using thermoelectric generators are wearable energy sources for ultra-low-power sensors, harvesters using small temperature gradients as energy sources, and on-chip thermoelectric harvesters [6][7][8][9][10][11][12][13][14][15][16].
There are two main approaches to the analysis and computer simulation of processes taking place in thermoelectric products. The first approach involves the use of finite element methods. This method allows multi-domain simulation with very high accuracy [12,14,17,18]. This method allows simulating an on-chip thermoelectric harvester, considering geometry and technology. At the same time, the method requires large computational resources and is mostly used to calculate a  Figure 1 depicts the three main types of thermoelectric elements. All the thermoelectric elements shown in the figure have a similar structure but differ in external dimensions, topology, geometry, and the number of elements in one module. The thermoelectric pair, indicated in the figure by the number 4 and the letters p and n, provides thermoelectric conversion. It is an integral part of the thermoelectric generator, consists of dissimilar materials, usually a pair of semiconductors of the p-and n-type with a pronounced Seebeck coefficient. About 30 of the most popular thermoelectric materials are described in [26]. The most common applications for each of the materials are also reported there. In the current work, materials that have properties close to those of polysilicon and (Bi 0.5 Sb 0.5 ) 2 Te 3 are used in the examples. Each of the "legs" may have a different geometry and may even be composed of segments from several different materials [27]. The thermoelectric "legs" are electrically connected and also connected to the other thermoelectric pairs in the module using conductive contacts, indicated in Figure 1 by the Number 3 and highlighted in red. A thermoelectric pair, together with the contact pads, are enclosed between the elements, Number 2, which provide electrical insulation, good thermal conductivity, and the mechanical rigidity of the entire structure. The thermal gradient over the thermoelectric pair leads to the appearance of electromotive force (EMF). The source of the thermal gradient is the difference between the temperature of the heat source and the ambient temperature T amb . In most cases, conductive heat transfer occurs between the thermoelectric module and the heat source, while convective heat transfer takes place between the thermoelectric module and T amb . through heatsink at Number 1.

Background on Thermoelectric Harvesters
Micromachines 2020, 11, x FOR PEER REVIEW 3 of 20 Figure 1 depicts the three main types of thermoelectric elements. All the thermoelectric elements shown in the figure have a similar structure but differ in external dimensions, topology, geometry, and the number of elements in one module. The thermoelectric pair, indicated in the figure by the number 4 and the letters p and n, provides thermoelectric conversion. It is an integral part of the thermoelectric generator, consists of dissimilar materials, usually a pair of semiconductors of the pand n-type with a pronounced Seebeck coefficient. About 30 of the most popular thermoelectric materials are described in [26]. The most common applications for each of the materials are also reported there. In the current work, materials that have properties close to those of polysilicon and (Bi0.5Sb0.5)2Te3 are used in the examples. Each of the "legs" may have a different geometry and may even be composed of segments from several different materials [27]. The thermoelectric "legs" are electrically connected and also connected to the other thermoelectric pairs in the module using conductive contacts, indicated in Figure 1 by the Number 3 and highlighted in red. A thermoelectric pair, together with the contact pads, are enclosed between the elements, Number 2, which provide electrical insulation, good thermal conductivity, and the mechanical rigidity of the entire structure. The thermal gradient over the thermoelectric pair leads to the appearance of electromotive force (EMF). The source of the thermal gradient is the difference between the temperature of the heat source and the ambient temperature Tamb. In most cases, conductive heat transfer occurs between the thermoelectric module and the heat source, while convective heat transfer takes place between the thermoelectric module and Tamb. through heatsink at Number 1. In a traditional thermoelectric element shown in Figure 1a, the size of the thermoelectric pair is dominant. Thus, the electrical and thermal resistances of the pair are dominant too. The electric contact resistance ((3) in Figure 1) and the thermal resistance of the insulating layers ((2) in Figure 1) are often assumed as negligible or less significant compared to the thermal and electrical properties of the pair itself [21,27,28]. A planar thermoelectric element shown in Figure 1b differs from the traditional one in that heat flows along its "legs" orthogonal to the direction of the general temperature gradient. With this arrangement, secondary temperature gradients along the "legs" can occur, which means that the temperature distribution along the "legs" becomes more complicated [29]. A thermoelectric pair in a planar on-chip thermoelectric element, shown in Figure 1c has a micrometric size [8][9][10]17,18]. Millions of thermoelectric pairs can be placed in one module. Such pairs are usually located on a silicon substrate. Special cavities in the substrate serve to thermally isolate the common contact between the p and n "legs". There may be air or a vacuum inside the cavities. The heat-conducting elements above and below the thermoelectric pair (Number 2 in Figure  1) transmute into thermal extenders, one of them is a silicon substrate and the other is a pillar consisting of metal layers and vias, whose dimensions significantly exceed the dimensions of the thermoelectric pair. This means that thermal processes outside the thermoelectric pair become no less significant than the processes within the pair itself. Besides, thermal gradients within the substrate In a traditional thermoelectric element shown in Figure 1a, the size of the thermoelectric pair is dominant. Thus, the electrical and thermal resistances of the pair are dominant too. The electric contact resistance ((3) in Figure 1) and the thermal resistance of the insulating layers ((2) in Figure 1) are often assumed as negligible or less significant compared to the thermal and electrical properties of the pair itself [21,27,28]. A planar thermoelectric element shown in Figure 1b differs from the traditional one in that heat flows along its "legs" orthogonal to the direction of the general temperature gradient. With this arrangement, secondary temperature gradients along the "legs" can occur, which means that the temperature distribution along the "legs" becomes more complicated [29]. A thermoelectric pair in a planar on-chip thermoelectric element, shown in Figure 1c has a micrometric size [8][9][10]17,18]. Millions of thermoelectric pairs can be placed in one module. Such pairs are usually located on a silicon substrate. Special cavities in the substrate serve to thermally isolate the common contact between the p and n "legs". There may be air or a vacuum inside the cavities. The heat-conducting elements above and below the thermoelectric pair (Number 2 in Figure 1) transmute into thermal extenders, one of them is a silicon substrate and the other is a pillar consisting of metal layers and vias, whose dimensions significantly exceed the dimensions of the thermoelectric pair. This means that Micromachines 2020, 11, 574 4 of 19 thermal processes outside the thermoelectric pair become no less significant than the processes within the pair itself. Besides, thermal gradients within the substrate can lead to different thermal conditions for different pairs. In some cases, this process must also be considered [17]. In addition to this, complementary metal-oxide-semiconductor (CMOS) or micro-electro-mechanical systems (MEMS) technology sometimes does not allow for the use of well-studied materials, requiring the trial of more recently developed ones. The thermal and electrical properties of such new materials, as well as their mutual influence, are still being researched. Many of the existing models of thermoelectric devices that reliably work with traditional thermoelectric elements do not take into account the asymmetry of heat distribution inherent in planar thermoelectric elements, the dependence of the thermal and electrical properties of materials on local temperature, and the inhomogeneity of temperature conditions for thermoelectric elements located at a considerable distance from each other. The spatial 3D model is required to analyze the operation of such thermoelectric modules.

Modeling
The spatial numerical solution of heat conduction in solids by the finite difference method includes dividing the body into several elemental volumes, compiling equations of temperature-heat flow at its boundaries and internal energy conversions. This method is described in detail in literature on heat transfer, for example [25]. An example of such an elemental volume is shown in Figure 2. The elemental volume is located in the Cartesian coordinate system. The faces of the elemental volume are named according to their location as x 1 , x 2 , y 1 , y 2 , z 1 , and z 2 correspondingly. The temperatures of the corresponding faces and the heat fluxes flowing in or out of the elemental volume have corresponding subscripts.

143
The law of energy conservation allows us to argue that the total heat entering the elemental 144 volume leads to a change in the internal energy of this elemental volume. From this statement, the 145 heat equation is derived, where is a value of the heat flowing into or out of the elemental volume through the corresponding 147 face, j means the coordinate axis, n means number 1 or 2 of the face, (see Figure. 2), is the heat 148 generated within the elementary volume, and E is the internal energy of the elemental volume.

149
The Fourier's Law of heat conduction along the x-axis has the following differential form: where Δ = 2 − 1 , and = , , . The is a temperature of the mass of the elemental volume.

154
For simplicity, it is assumed that the mass of the elemental volume is concentrated in the middle of 155 the volume. The law of energy conservation allows us to argue that the total heat entering the elemental volume leads to a change in the internal energy of this elemental volume. From this statement, the heat equation is derived, where q ji is a value of the heat flowing into or out of the elemental volume through the corresponding face, j means the coordinate axis, n means number 1 or 2 of the face, (see Figure 2), q g is the heat generated within the elemental volume, and E is the internal energy of the elemental volume. The Fourier's Law of heat conduction along the x-axis has the following differential form: where k is thermal conductivity, [W·m −1 ·K −1 ], A x is the area of the face x 1 . Rewriting (2) in the form of differences for corresponding axes, and expressing internal energy of the elemental volume as the specific heat of material c p , [J·kg −1 · K −1 ] times its temperature, we get: where ∆T j = T j2 − T j1 , and j = x, y, or z. The T C is a temperature of the mass of the elemental volume. For simplicity, it is assumed that the mass of the elemental volume is concentrated in the middle of the volume.

Electrical Analogy
In this article, we will first use the equivalent circuit developed in [24], then expand it to include an electrical domain for mutual simulation of thermal and electrical processes, and finally add thermoelectric process to the expanded circuit. The thermoelectric analogies used in [19,21,24,25] are listed in Table 1.
k · ∆z ∆y·∆x : Thermal resistances of the elemental volume in x, y and z directions, [K/W] q j = ∆T j /Θ j : Fourier's law along axis j i j = ∆ϕ j R j : Ohm's law along axis j q C = C T · dT C dt : change of internal energy i C = C· dV C dt : charge of a capacitor For the sake of reconstructing the equivalent circuit provided in [24], let us rewrite Equation (3) using the relations listed in Table 1.
An equivalent circuit corresponding to (4) is shown in Figure 3, thermal domain. A mass of the elemental volume is presented as an equivalent capacitor C T . One terminal of the capacitor is connected to the central junction, which is separated from each of the faces by the corresponding thermal resistance where R j1 = R j2 = 0.5 R j , where subscript j means x−, y− or z−axis name. The second terminal of the capacitor is connected to the common point of the thermal domain (0, kelvin). An equivalent current source q g that represents the heat generated within the elemental volume is also connected between zero potential and a central junction. An initial condition (IC) of the equivalent capacitor is the initial temperature of the elemental volume. The boundary conditions may be set by the current source as an equivalent of heat flow, as an equivalent of temperature, by a voltage source, or by contact with neighboring elemental volume. Micromachines 2020, 11, x FOR PEER REVIEW 6 of 20 If the heat qg generated inside the elemental volume is the result of the Joule heating, then this process can also be included in the model for simulation by adding an electrical domain to the model: see Figure 3, electrical domain. The electrical model includes only six resistors, equivalent to the electrical resistance of the elemental volume in each dimension. 1 = 2 = 2 , see Table 1. The value of the Joule heating is equal to a voltage drop over each resistor ( ) times current ( ) through this resistor.

Extending the Model by Introducing the Thermoelectric Processes
The model shown in Figure 3 should be expanded to consider the Seebeck and Peltier thermoelectric effects. In an expanded form, this model will be suitable for modeling thermoelectric elements.
The essence of the Seebeck thermoelectric process is that in the presence of a thermal gradient over the conductor, an EMF appears in the conductor. The magnitude of the EMF is proportional to the temperature drop. This phenomenon is explained by the fact that the greater the energy possessed by the free charge carriers in a conductive material, the higher the temperature of their location zone. As a result, charge carriers drift from the warmer zone of the conductive material towards the colder zone due to diffusion. The coefficient that relates the EMF to the temperature gradient is called the Seebeck coefficient ( , V • K −1 ). In some materials, the Seebeck coefficient is significantly higher than in others. The Seebeck coefficient can be both negative in the case when the charge carriers are negative (electrons) and positive if the charge carriers have a positive sign (holes). A thermoelectric pair, such as the one in Figure 1, consists of two dissimilar materials with Seebeck coefficients of different signs (p-and n-type semiconductors). In the presence of the same temperature gradient, such materials generate EMFs with the opposite sign. These EMFs are summed up due to the galvanic contact between the materials, as shown in Figure 4.
Since the Seebeck effect is a property of the material, it can be described for an elemental volume. Assume that a temperature gradient is applied to the elemental volume along the global axis J. The temperatures at the boundaries of the volume are 1 and 2 . In this case, the EMF arising at the electrical terminals of this volume is equal to: If the heat q g generated inside the elemental volume is the result of the Joule heating, then this process can also be included in the model for simulation by adding an electrical domain to the model: see Figure 3, electrical domain. The electrical model includes only six resistors, equivalent to the electrical resistance of the elemental volume in each dimension. R j1 = R j2 = R j 2 , see Table 1. The value of the Joule heating is equal to a voltage drop over each resistor (v R jn ) times current (i jn ) through this resistor.

Extending the Model by Introducing the Thermoelectric Processes
The model shown in Figure 3 should be expanded to consider the Seebeck and Peltier thermoelectric effects. In an expanded form, this model will be suitable for modeling thermoelectric elements.
The essence of the Seebeck thermoelectric process is that in the presence of a thermal gradient over the conductor, an EMF appears in the conductor. The magnitude of the EMF is proportional to the temperature drop. This phenomenon is explained by the fact that the greater the energy possessed by the free charge carriers in a conductive material, the higher the temperature of their location zone. As a result, charge carriers drift from the warmer zone of the conductive material towards the colder zone due to diffusion. The coefficient that relates the EMF to the temperature gradient is called the Seebeck coefficient (α, V·K −1 ). In some materials, the Seebeck coefficient is significantly higher than in others. The Seebeck coefficient can be both negative in the case when the charge carriers are negative (electrons) and positive if the charge carriers have a positive sign (holes). A thermoelectric pair, such as the one in Figure 1, consists of two dissimilar materials with Seebeck coefficients of different signs (p-and n-type semiconductors). In the presence of the same temperature gradient, such materials generate EMFs with the opposite sign. These EMFs are summed up due to the galvanic contact between the materials, as shown in Figure 4.
Since the Seebeck effect is a property of the material, it can be described for an elemental volume. Assume that a temperature gradient is applied to the elemental volume along the global axis J. The temperatures at the boundaries of the volume are T J1 and T J2 . In this case, the EMF arising at the electrical terminals of this volume is equal to: Thus, we can say that the electric potential φ j1 appears on the face j 1 of the elemental volume under the influence of the temperature T j1 , and the electric potential φ j2 appears on the opposite face j 2 under the influence of the temperature T j2 . A pair of oppositely directed behavioral voltage sources v s1 and v s2 , equal in magnitude to the thermoelectric potentials φ j1 and φ j2 , respectively, are shown in the electrical region of Figure 4. If T j1 > T j2 , then the electric potential on φ j1 surface j 1 is higher than φ j2 on surface j 2 ; thus, EMF is nonzero.
Micromachines 2020, 11, x FOR PEER REVIEW 7 of 20 Thus, we can say that the electric potential 1 appears on the face j1 of the elemental volume under the influence of the temperature 1 , and the electric potential 2 appears on the opposite face is the ohmic resistance of the elemental volume on j-axis and is its thermal resistance on j-axis.
When a load is connected to electrical terminals, an electric current ij flows through the volume. In the circuit on Figure 4, the electric current passes through the voltage source 2 in a positive direction, and through the source 2 in а negative direction. In other words, a certain amount of electric power equal to the voltage drop over the voltage source 2 multiplied by the current ij is converted (pumped) from electrical power at the electrical domain into heat in the thermal domain in face j2. The electrical negative power is raised at voltage source 1 when the current ij flows through it in the negative direction. This power is pumped from the thermal domain at face j1 to the electrical domain. The thermoelectric effect, which leads to the absorption of heat on one surface of the material and the emission of heat on the other surface as a result of an electric current in the material, is called the Peltier effect. Formally, the Peltier effect is expressed as follows: where is Number 1 or 2 of the face of elemental volume along the corresponding axis.
Using the same approach, to emulate Peltier and Seebeck thermoelectric effects, we can add the corresponding elements into the spatial equivalent circuit shown in Figure 3. Figure 5 depicts a modified spatial equivalent circuit of an elemental volume suitable for simulating thermoelectric processes.  Figure 3. The global axis j is used instead of x-, yor z-axes. R j is the ohmic resistance of the elemental volume on j-axis and Θ j is its thermal resistance on j-axis.
When a load is connected to electrical terminals, an electric current i j flows through the volume. In the circuit on Figure 4, the electric current passes through the voltage source v s2 in a positive direction, and through the source v s2 in a negative direction. In other words, a certain amount of electric power equal to the voltage drop over the voltage source v s2 multiplied by the current i j is converted (pumped) from electrical power at the electrical domain into heat in the thermal domain in face j 2 . The electrical negative power is raised at voltage source v s1 when the current i j flows through it in the negative direction. This power is pumped from the thermal domain at face j 1 to the electrical domain. The thermoelectric effect, which leads to the absorption of heat on one surface of the material and the emission of heat on the other surface as a result of an electric current in the material, is called the Peltier effect. Formally, the Peltier effect is expressed as follows: where n is Number 1 or 2 of the face of elemental volume along the corresponding axis.
Using the same approach, to emulate Peltier and Seebeck thermoelectric effects, we can add the corresponding elements into the spatial equivalent circuit shown in Figure 3. Figure 5 depicts a modified spatial equivalent circuit of an elemental volume suitable for simulating thermoelectric processes.

Temperature Coefficient
The two-domain equivalent circuit of the elemental volume of the material shown in Figure 5 already allows simulations of thermoelectric processes in materials. However, this circuit has one significant drawback: it employs physical parameters, such as electrical resistivity, thermal conductivity, and Seebeck coefficient, as constant values. In practice, all of these physical quantities are subject to variations due to temperature. These changes are not negligible. An increase in the electrical resistance of a material with temperature leads to significant losses in output power, and the temperature dependence of the Seebeck coefficient on temperature significantly affects the output voltage of a thermoelectric generator. This phenomenon is known as the Thomson thermoelectric effect. This effect can be neglected in the analysis of a traditional thermoelectric generator, but with planar or on-chip generators, this effect should not be neglected.

Temperature Coefficient
The two-domain equivalent circuit of the elemental volume of the material shown in Figure 5 already allows simulations of thermoelectric processes in materials. However, this circuit has one significant drawback: it employs physical parameters, such as electrical resistivity, thermal conductivity, and Seebeck coefficient, as constant values. In practice, all of these physical quantities are subject to variations due to temperature. These changes are not negligible. An increase in the electrical resistance of a material with temperature leads to significant losses in output power, and the temperature dependence of the Seebeck coefficient on temperature significantly affects the output voltage of a thermoelectric generator. This phenomenon is known as the Thomson thermoelectric effect. This effect can be neglected in the analysis of a traditional thermoelectric generator, but with planar or on-chip generators, this effect should not be neglected.
As a simplified function of the dependence of a certain parameter X(T) on temperature, one can use the first-order function, expressed using the thermal coefficient: where * is a value of parameter X(T) measured at reference temperature * , T is an actual temperature, and x is a temperature coefficient. Thus, we can express the necessary parameters of materials in Form (8) The physical parameters of most of the materials, as well as corresponding thermal coefficients, are available in the literature.
Temperature is one of the variables to be simulated: we want to define a parameter of a model as a function of temperature. The values of the electrical resistors and the thermal resistors represented in the model by the equivalent electrical resistors used in the model are temperature invariable. The temperature dependent resistors can be modeled using behavioral sources.
Ohm's law states the following: the voltage drop across the resistor is proportional to the current: As a simplified function of the dependence of a certain parameter X(T) on temperature, one can use the first-order function, expressed using the thermal coefficient: where X * is a value of parameter X(T) measured at reference temperature T * , T is an actual temperature, and x is a temperature coefficient. Thus, we can express the necessary parameters of materials in Form (8) ρ The physical parameters of most of the materials, as well as corresponding thermal coefficients, are available in the literature.
Temperature is one of the variables to be simulated: we want to define a parameter of a model as a function of temperature. The values of the electrical resistors and the thermal resistors represented in the model by the equivalent electrical resistors used in the model are temperature invariable. The temperature dependent resistors can be modeled using behavioral sources.
Ohm's law states the following: the voltage drop across the resistor is proportional to the current: Thus, the voltage drop across the resistor due to current can be emulated with a current-controlled voltage source (CCVS). When using a CCVS, the output port of the voltage source is connected to the branch of the circuit instead of the resistor. The voltage drop, in this case, is set to be equal to the product of current flowing through the branch and the resistance value. The resistance value can be specified as a number (temperature independent value) as well as an Equations (9) and (10).
The temperature dependence of the Seebeck coefficient can be introduced into the model by replacing the constant value of the α in the v jn sources with Equation (11).
An example of a hierarchical block representing the equivalent circuit of an elemental volume from Figure 5 using the described technique is shown in Figure 6. The block has 12 bidirectional connectors named as f i jn for the electrical domain and t jn for the thermal domain. Subscript j is a corresponding x, y, or z axis name, and n is the surface number.
The temperature dependence of the Seebeck coefficient can be introduced into the model by replacing the constant value of the in the sources with Equation (11) An example of a hierarchical block representing the equivalent circuit of an elemental volume from Figure 5 using the described technique is shown in Figure 6. The block has 12 bidirectional connectors named as for the electrical domain and for the thermal domain. Subscript is a corresponding , , axis name, and is the surface number. Figure 6. An equivalent circuit for simulation of thermoelectric processes in elemental volume. The equivalent circuit is designed as a hierarchical block to simulate in LTSpice software [30]. Parameters , ℎ, ℎ, ℎ, 1, 0, 1, 1, ℎ, ℎ, ℎ, and means , In Figure 7 an example of the schematic symbol proposed for such a type of equivalent circuit is shown. The symbol of the hierarchical block can be of any shape. We propose to make it look like a cube, in order to emphasize the association of the equivalent circuit with the 3D element of the spatial structure. Each face of the cube has two terminals: one terminal ( ) belongs to the electrical domain and the second terminal ( ) belongs to the thermal domain.  Figure 6 for LTSpice [30]. The equivalent circuit is designed as a hierarchical block to simulate in LTSpice software [30]. Parameters ct, Rxh, Ryh, Rzh, r1, se0, se1, k1, THxh, THyh, THzh, and tre f means C T , R x /2, R y /2, R z /2, ρ Ω0 , α * , α 0 , k 0 , Θ x /2, Θ y /2, Θ z /2, and T * , respectively. The resistors R 1 − R 12 in the scheme have insignificant values and serve to prevent the problem of convergence during circuit simulation. In Figure 7 an example of the schematic symbol proposed for such a type of equivalent circuit is shown. The symbol of the hierarchical block can be of any shape. We propose to make it look like a cube, in order to emphasize the association of the equivalent circuit with the 3D element of the spatial structure. Each face of the cube has two terminals: one terminal ( f i jn ) belongs to the electrical domain and the second terminal (t jn ) belongs to the thermal domain.
The temperature dependence of the Seebeck coefficient can be introduced into the model by replacing the constant value of the in the sources with Equation (11) An example of a hierarchical block representing the equivalent circuit of an elemental volume from Figure 5 using the described technique is shown in Figure 6. The block has 12 bidirectional connectors named as for the electrical domain and for the thermal domain. Subscript is a corresponding , , axis name, and is the surface number. Figure 6. An equivalent circuit for simulation of thermoelectric processes in elemental volume. The equivalent circuit is designed as a hierarchical block to simulate in LTSpice software [30]. Parameters , ℎ, ℎ, ℎ, 1, 0, 1, 1, ℎ, ℎ, ℎ, and means , In Figure 7 an example of the schematic symbol proposed for such a type of equivalent circuit is shown. The symbol of the hierarchical block can be of any shape. We propose to make it look like a cube, in order to emphasize the association of the equivalent circuit with the 3D element of the spatial structure. Each face of the cube has two terminals: one terminal ( ) belongs to the electrical domain and the second terminal ( ) belongs to the thermal domain.  Figure 6 for LTSpice [30].  Figure 6 for LTSpice [30].
The electrical domain of one block may be joined to the electrical domain of other blocks or electrical elements such as an electrical load or a current or voltage source. The thermal domain of one block may be connected to the thermal domain of another block, or equivalent electrical elements representing a temperature source, heat sources, or thermal impedances. The thermal and electrical domains should remain isolated from each other.

Examples of Simulations
In this part, examples of heat-conducting simulation using the proposed method for solving various spatial problems are demonstrated. These examples relate to solving one-, two-, and three-dimensional heat conduction problems in solids, to simulating thermoelectric processes in a single material and in a thermoelectric pair. Further, this section analyzes the behavior of a group of thermoelectric elements under similar conditions, macro models including millions of thermoelectric pairs located on a single substrate, and peripheral devices such as a radiator and an external electrical load.
The materials used in the examples have physical properties close to those of the materials normally used for building on-chip thermoelectric generators in CMOS technology. Some properties of the materials are deliberately exaggerated in order to visualize certain physical phenomena, such as the Thomson effect. The topology of the sample thermoelectric generator is also demonstrated in a simplified form with approximate geometric dimensions.

Heat Conduction Simulation
As the first example, consider the heat conduction problem shown in Figure 8. The plate is made of a material whose properties are close to those of aluminum. The plate is thermally insulated from above, from below, and partially from the front (insulation is shown in white in the picture). Prior to the simulation, the plate had an initial temperature equal to the ambient temperature T amb . At the beginning of the simulation, boundary conditions were applied to the plate in the form of temperatures T 1 on the left face, temperature T 2 in the middle third of the front face, and temperature T 3 on the right and back faces. At the beginning of the experiment, the front and back faces of the plate were electrically grounded. The electrical potential V 1 was supplied to the rear face of the plate with a delay t 1 relative to the start of the experiment.

Examples of Simulations
In this part, examples of heat-conducting simulation using the proposed method for solving various spatial problems are demonstrated. These examples relate to solving one-, two-, and threedimensional heat conduction problems in solids, to simulating thermoelectric processes in a single material and in a thermoelectric pair. Further, this section analyzes the behavior of a group of thermoelectric elements under similar conditions, macro models including millions of thermoelectric pairs located on a single substrate, and peripheral devices such as a radiator and an external electrical load.
The materials used in the examples have physical properties close to those of the materials normally used for building on-chip thermoelectric generators in CMOS technology. Some properties of the materials are deliberately exaggerated in order to visualize certain physical phenomena, such as the Thomson effect. The topology of the sample thermoelectric generator is also demonstrated in a simplified form with approximate geometric dimensions.

Heat Conduction Simulation
As the first example, consider the heat conduction problem shown in Figure 8. The plate is made of a material whose properties are close to those of aluminum. The plate is thermally insulated from above, from below, and partially from the front (insulation is shown in white in the picture). Prior to the simulation, the plate had an initial temperature equal to the ambient temperature Tamb. At the beginning of the simulation, boundary conditions were applied to the plate in the form of temperatures T1 on the left face, temperature T2 in the middle third of the front face, and temperature T3 on the right and back faces. At the beginning of the experiment, the front and back faces of the plate were electrically grounded. The electrical potential V1 was supplied to the rear face of the plate with a delay t1 relative to the start of the experiment. The parameters used in the simulation are collected in Table 2.  The parameters used in the simulation are collected in Table 2. Table 2. Parameters used in the simulation.

Type of Parameter Parameter Value Units
Dimensions Parameters of simulation Ambient air temperature, T amb 25 + 273.15 K T 1 T amb + 10 K T 2 T amb + 5 K T 3 T amb K V 1 delayed by 0.5 s 40 × 10 −3 V Figure 9 shows the electrical equivalent circuit of the problem depicted in Figure 8, performed in the schematic editor of the LTSpice circuits simulator software. The plate is divided into a mesh of 18 blocks, 6 along the x-axis and 3 along the y-axis, and one along the z-axis as shown in Figure 6. Parameters are listed out on the working page of the circuit editor. Temperatures are given as equivalent voltage sources.  Figure 9 shows the electrical equivalent circuit of the problem depicted in Figure 8, performed in the schematic editor of the LTSpice circuits simulator software. The plate is divided into a mesh of 18 blocks, 6 along the x-axis and 3 along the y-axis, and one along the z-axis as shown in Figure 6. Parameters are listed out on the working page of the circuit editor. Temperatures are given as equivalent voltage sources. The temperature response of each elemental volume of the plate to the boundary conditions and the internal heat generated due to electrical current that turns on after 0.5 delay is shown in Figure  10 in the form of curves of the equivalent voltage vs. time. The 1 V of the equivalent voltage corresponds to 1 °C. Figure 9. The electrical equivalent circuit in LTSpice schematic editor for the heat conduction problem from Figure 8. The circuit includes 18 blocks shown in Figure 6. The blocks are designated as X i j , where i is the row number and the j is a block's number in the row.
The temperature response of each elemental volume of the plate to the boundary conditions and the internal heat generated due to electrical current that turns on after 0.5 s delay is shown in Figure 10

Thermoelectric Couple
This section discusses an example of simulating the behavior of a planar on-chip thermoelectric couple Figure 11. Approximate dimensions have been selected based on [7,8,17,23,31]. The thermoelectric couple consists of a pair of semiconductor materials of p-and n-types (p-leg and n-leg

Thermoelectric Couple
This section discusses an example of simulating the behavior of a planar on-chip thermoelectric couple Figure 11.
Approximate dimensions have been selected based on [7,8,17,23,31]. The thermoelectric couple consists of a pair of semiconductor materials of p-and n-types (p-leg and n-leg in the figure). The couple is mounted on a silicon substrate with the cavity etched under the cold junction for thermal isolation between the cold and the hot edges of the thermoelectric legs. The hot external edges of the thermoelectric pair have thermal contact with the substrate (substrate 1 and substrate 2) and the substrate, in turn, is connected to a hot temperature source (T 1 ) at the bottom side of the construction. The cold internal edges of the thermoelectric legs are connected thermally to an upper metal layer (metal 2) through the heat-conducting extender made from metal layers and metallized vias. The extender is assumed for the sake of simplicity as a uniform material (metal 1). The upper metal layer is connected thermally to a cold temperature source (T 2 ). The properties of the materials used in the simulation are tabulated in Table 3. The materials used to model the thermoelectric pair are taken hypothetically because the model is generic. The properties of the materials are close to those of real materials, but we can vary them slightly to obtain more visual simulation results. For simplicity, suppose that there is a vacuum in the cavities and the thermal insulation is ideal. Otherwise, the properties of the insulating material should be added, as shown in [17].
The properties of the materials used in the simulation are tabulated in Table 3. Figure 12 shows the load curves (V-I and P-I) for the single thermoelectric pair being tested. The result was obtained using DC-sweep analysis for various values of the thermal coefficient α 0 of the Seebeck coefficient α for both p-and n-materials. This experiment demonstrates the Thomson effect on the output characteristics of the module. We remind that α 0 values are taken to demonstrate the effect only and are not physical properties of specific materials. Figure 12 shows the voltage-current and power-current characteristic curves of the single thermoelectric pair, shown in Figure 11. The power yield of one pair is small, but the idea of a thermoelectric harvester on a chip involves the use of many (up to several million) such pairs. All the pairs are normally divided into sub-units of several thousand pairs each, electrically connected in series. The groups can be interconnected in series or parallel, depending on the power requirements, [23]. In [17] it is demonstrated that thermoelectric pairs located at a large physical distance from each other are in slightly different temperature conditions. Nevertheless, the temperature conditions for pairs in close proximity to one another can be considered identical. Table 3. Properties of the materials used in the thermoelectric couple under simulation.

Metal 2
Dimensions L x = 141 µm, L y = 4 µm, and L z = 2 µm Meshing along x, y, and z axes 5 × 1 × 1 Thermodynamic properties see metal 1 Substrate 1 Dimensions L x = 4 µm, L y = 8 µm, and L z = 4 µm Meshing along x, y, and z axes 1 × 1 × 1 Thermodynamic properties  Figure 11 for various values of thermal coefficient of Seebeck constants for p-and n-types of thermoelectric materials. Figure 12 shows the voltage-current and power-current characteristic curves of the single thermoelectric pair, shown in Figure 11. The power yield of one pair is small, but the idea of a thermoelectric harvester on a chip involves the use of many (up to several million) such pairs. All the pairs are normally divided into sub-units of several thousand pairs each, electrically connected in series. The groups can be interconnected in series or parallel, depending on the power requirements, [23]. In [17] it is demonstrated that thermoelectric pairs located at a large physical distance from each other are in slightly different temperature conditions. Nevertheless, the temperature conditions for pairs in close proximity to one another can be considered identical.

Sub-Units' Simulation
As an example of the proposed method's usage for multiple pairs and sub-units, we analyze the operation of a thermoelectric harvester on a chip that includes 9 sub-units, one thousand pairs each. The elements of one sub-unit are considered to be located close enough to each other so that it can be said that their functioning is identical. In this case, there is no need to multiply the number of blocks. The model of a single pair simulates the behavior of all the pairs in the entire sub-unit.  Figure 11 for various values of thermal coefficient of Seebeck constants for p-and n-types of thermoelectric materials.
As an example of the proposed method's usage for multiple pairs and sub-units, we analyze the operation of a thermoelectric harvester on a chip that includes 9 sub-units, one thousand pairs each. The elements of one sub-unit are considered to be located close enough to each other so that it can be said that their functioning is identical. In this case, there is no need to multiply the number of blocks. The model of a single pair simulates the behavior of all the pairs in the entire sub-unit.
Electrically, all couples in a sub-unit are connected in series. This means that their total output voltage is equal to the output voltage of a single couple times the number of couples (noc) in the sub-unit. The common electric current of all couples of the sub-unit is the same for all pairs and depends on the load. The electrical domain in Figure 13a demonstrates the implementation of these principles using a pair of behavioral sources of current and voltage. A controlled current source provides an output current of the couple equal to the load current and a controlled voltage source provides the load voltage equal to the voltage of a single couple times the number of couples (noc) of the sub-unit.
The elements of one sub-unit are considered to be located close enough to each other so that it can be said that their functioning is identical. In this case, there is no need to multiply the number of blocks. The model of a single pair simulates the behavior of all the pairs in the entire sub-unit.
Electrically, all couples in a sub-unit are connected in series. This means that their total output voltage is equal to the output voltage of a single couple times the number of couples ( ) in the subunit. The common electric current of all couples of the sub-unit is the same for all pairs and depends on the load. The electrical domain in Figure 13a demonstrates the implementation of these principles using a pair of behavioral sources of current and voltage. A controlled current source provides an output current of the couple equal to the load current and a controlled voltage source provides the load voltage equal to the voltage of a single couple times the number of couples ( ) of the sub-unit.
(a) (b) Figure 13. Simulation of a 1k-couples sub-unit using the model of a single couple shown in Figure 11.  The parallel connection of the couples in a thermal domain can be imitated similarly, using two pairs of interdependent sources, see Figure 13a. In this case, it should be remembered that the amount of heat flowing into the couple from the hot side is not equal to the amount of heat flowing from the cold side, since part of the heat is converted to electricity and part spent on heat loss. Thus, the thermal terminals h and c should be considered independent ports, where the common terminal is 0 kelvin. For each of these thermal ports, the corresponding behavioral voltage source sets the equivalent temperature equal to the temperature of the source, and the behavioral current source multiplies the equivalent heat flux of the temperature source by the number of couples (NOC). The more pairs working in parallel, the higher the heat flow through the temperature sources.

Macro-Level Simulation
Of particular interest are joint simulations of the whole chip, which include several subunits that receive heat from sources through heat exchangers and thermal interfaces, and loaded with a converter that converts the output voltage of the harvester to one of the standard voltages. To avoid overloading this article with multiple, even though interesting simulations, we offer the simplified macro simulation that will demonstrate the capabilities of the method. If desired, this simulation can be complicated by simply adding elements.
A fragment of a harvester on a chip is shown in Figure 14a. One of the nine subunits gets heat from the heat source and releases it into the surrounding air. Convection heat transfer is carried out using a heatsink. The thermal resistance between the substrate and the heat source is much less than the thermal resistance of the heatsink, so it can be neglected. A review of heatsinks shows that a heatsink with a base of 25 mm × 25 mm demonstrates thermal resistance of about 30 K·W −1 at zero airflow speed [32]. This means that the fragment of the heatsink (1/9 of its area) has approximately 9 times greater thermal resistance. Figure 14b shows the equivalent scheme in LTSpice and its load characteristics are shown in Figure 14c. using a heatsink. The thermal resistance between the substrate and the heat source is much less than the thermal resistance of the heatsink, so it can be neglected. A review of heatsinks shows that a heatsink with a base of 25 mm × 25 mm demonstrates thermal resistance of about 30 K • W −1 at zero airflow speed [32]. This means that the fragment of the heatsink (1/9 of its area) has approximately 9 times greater thermal resistance. Figure 14b shows the equivalent scheme in LTSpice and its load characteristics are shown in Figure 14c.

Validation of the Method and Comparison with a Proved Methodology
The method of spatial equivalent circuit modeling proposed in this study is derived from the finite difference method, in which the differential equations for an elemental volume of a threedimensional grid are expressed in the form of equivalent electrical circuits. A similar technique is often used to solve thermal conductivity problems [21,24,25]. The well-known equivalent circuit widely used to solve thermal conductivity problems is expanded in this study. The equations of thermoelectric processes are added to the Fourier equation of the elemental volume and are also expressed as elements of an equivalent electrical circuit.

Validation of the Method and Comparison with a Proved Methodology
The method of spatial equivalent circuit modeling proposed in this study is derived from the finite difference method, in which the differential equations for an elemental volume of a three-dimensional grid are expressed in the form of equivalent electrical circuits. A similar technique is often used to solve thermal conductivity problems [21,24,25]. The well-known equivalent circuit widely used to solve thermal conductivity problems is expanded in this study. The equations of thermoelectric processes are added to the Fourier equation of the elemental volume and are also expressed as elements of an equivalent electrical circuit.
The objective of this section is to validate the legitimacy of using the proposed method in its expanded form. For this purpose, we chose a paper "Multiphysics simulation of thermoelectric systems-modeling of Peltier-cooling and thermoelectric generation" [33], in which the model of thermoelectric Peltier cooling and thermoelectric generation processes was simulated by means of the finite element method using COMSOL Multiphysics software (3.3a, COMSOL, Inc., Burlington, MA, USA). Simulations for steady-state processes under different boundary conditions and electric current values, as well as simulations of transients provided in the article, were reconstructed using the proposed method.
A simple cooler geometry consisted of one p-type (Bi 0.5 Sb 0.5 ) 2 Te 3 semiconductor element 1 mm × 1 mm × 5.8 mm in size (see Figure 15a). It was contacted by two copper electrodes 0.1 mm in thickness. The thermal, electrical, and thermoelectric properties of the material as function of temperature are tabulated in Table 4. The properties of copper are assumed temperature independent. The material properties of the copper, from [33], were as follows: α c = 6.5 [µV·K −1 ], σ c = 5.9 × 10 8 , [S·m −1 ], and λ c = 350, [W·m −1 ·K −1 ].

Example 1: Steady-State Operation
The boundary conditions were set to 0 V and 0 • C at the base of the lower electrode. Adiabatic boundary conditions were taken on all other surfaces. At the top of the upper electrode, a current of 0.7 A was applied. The top side of the element under test was thermally isolated (thermally unloaded) or thermally loaded by 10 mW thermal load in rotation. The results of the simulation using the finite element method in COMSOL Multiphysics software presented in [33] are shown in the picture by solid gray lines. The graphics from [33] were digitized using Engauge digitizer software [34]. Figure 15b demonstrates the cold side temperatures versus current without and with 10 mW heat load, respectively.  Figure 16 demonstrates the reaction of the thermoelectric system to the current pulse in the time domain. The thermoelectric element described above is operated at 0.6 ampere with no heat load at steady-state temperatures. One second after the beginning of the calculation, the current is increased to 0.8 A linearly within 0.1 s. After 4 seconds the amperage decreases back to 0.6 A linearly within 0.1s (see the upper curve in Figure 16). The solid gray line shows the transient cold side temperature of the element, with temperature-dependent material parameters simulated in COMSOL Multiphysics from [33]. The dashed line shows the results of transient simulation of the element by means of LTSpice circuit simulation software using the method proposed in the current study. A slight deviation of the obtained results from the results that we accepted as a reference is explained by the fact that the mesh we used in the equivalent circuit simulation has only 12 blocks. The number of nodes used to simulate finite element models in specialized software is usually significantly higher. Nevertheless, the results of calculations performed by two different methods have very close values and show the nature of the temperature change as a reaction to excitation in a similar way.

Discussion and Conclusions
The technique described above allows the simulation of thermal and thermoelectric processes together with electrical ones. The method is based on the creation of a two-domain (electrical and  In order to simulate the system using the method proposed in this study, the twelve hierarchical blocks of the elemental volumes that are shown in Figure 6 were used. (Ten blocks for the thermoelectric element model and two for the copper contacts). The dependence of the properties of the material on temperature was performed using the lookup tables in accordance with the values of Table 4. The blocks were connected in series along the z-axis both electrically and thermally. The electrical domain of the first block was connected to a 0 V electrical potential, and the thermal domain of the first block was connected to a source of the equivalent temperature of 273.15 K (0 • C). The electrical domain of the last block was connected to an electric current source and the thermal domain of the last block was loaded by the equivalent heat source.
Simulation results are shown by dashed lines in Figure 15b,c together with results obtained in [33]. Comparing the results, we can conclude that in the steady state, the proposed method gives practically identical results to the finite-element method. Figure 16 demonstrates the reaction of the thermoelectric system to the current pulse in the time domain. The thermoelectric element described above is operated at 0.6 ampere with no heat load at steady-state temperatures. One second after the beginning of the calculation, the current is increased to 0.8 A linearly within 0.1 s. After 4 seconds the amperage decreases back to 0.6 A linearly within 0.1 s (see the upper curve in Figure 16). The solid gray line shows the transient cold side temperature of the element, with temperature-dependent material parameters simulated in COMSOL Multiphysics from [33]. The dashed line shows the results of transient simulation of the element by means of LTSpice circuit simulation software using the method proposed in the current study. A slight deviation of the obtained results from the results that we accepted as a reference is explained by the fact that the mesh we used in the equivalent circuit simulation has only 12 blocks. The number of nodes used to simulate finite element models in specialized software is usually significantly higher. Nevertheless, the results of calculations performed by two different methods have very close values and show the nature of the temperature change as a reaction to excitation in a similar way.

Example 2. Transient Operation
domain. The thermoelectric element described above is operated at 0.6 ampere with no heat load at steady-state temperatures. One second after the beginning of the calculation, the current is increased to 0.8 A linearly within 0.1 s. After 4 seconds the amperage decreases back to 0.6 A linearly within 0.1s (see the upper curve in Figure 16). The solid gray line shows the transient cold side temperature of the element, with temperature-dependent material parameters simulated in COMSOL Multiphysics from [33]. The dashed line shows the results of transient simulation of the element by means of LTSpice circuit simulation software using the method proposed in the current study. A slight deviation of the obtained results from the results that we accepted as a reference is explained by the fact that the mesh we used in the equivalent circuit simulation has only 12 blocks. The number of nodes used to simulate finite element models in specialized software is usually significantly higher. Nevertheless, the results of calculations performed by two different methods have very close values and show the nature of the temperature change as a reaction to excitation in a similar way. Figure 16. Transient calculation of Peltier-supercooling. A short current pulse (upper magenta line) generates temporarily lower cold side temperatures. The gray solid line demonstrates results of the finite element simulation from [33], the dashed line corresponds to a simulation by the proposed method.

Discussion and Conclusions
The technique described above allows the simulation of thermal and thermoelectric processes together with electrical ones. The method is based on the creation of a two-domain (electrical and

Discussion and Conclusions
The technique described above allows the simulation of thermal and thermoelectric processes together with electrical ones. The method is based on the creation of a two-domain (electrical and thermal) equivalent circuit of elemental volume with energy interrelations between the domains. The equations in both domains are solved for each of the three spatial axes; therefore, this method allows the simulation of processes for which the spatial arrangement matters.
The article is illustrated with examples of using the proposed method: 1. Modeling for three-dimensional time-domain analysis of the heat distribution in solids due to boundary conditions variation and internal heat generation. 2. Steady-state analysis of a planar on-chip thermoelectric couple, and a sub-unit including 1k couples. 3. Analysis of a thermoelectric sub-unit connected to heat/cold source through the heatsink. 4. The influence of the dependence of the thermodynamic and electrical parameters of the material on the local temperature in each of the mesh nodes is demonstrated in the example of the Thomson effect when the thermoelectric coefficient changes its value from node to node of the mesh, depending on the temperature of each of the nodes.
A hierarchical block, which includes a two-domain equivalent elemental volume scheme, is the basic element for creating models. Thermodynamic and electrical parameters of the material are the parameters of the corresponding hierarchical block. The more blocks involved in the simulation, the higher the spatial resolution, and therefore the accuracy.
Any software for electrical circuits simulation can provide the simulations using the proposed method. In this work, LTSpice, the SPICE-based simulator, was used to produce examples. The method allows for the use of all the main types of simulations used in engineering: steady-state (bias point, DC-sweep, etc.) analysis, time-domain (transient) analysis, as well as frequency domain (AC) analysis. Any external electrical components (such as load, PWM converter, energy storage, close-loop voltage regulator, etc.) can be connected to the ports of the electrical domain of the equivalent circuit and participate in the joint simulation.
A comparison of the results of simulations performed using the proposed method with similar simulations performed by other researchers using the finite element method shows a perfect match in a large number of examples and a relatively small difference in the rest of the examples.
One drawback of the method is the cumbersomeness of equivalent circuits, entailing high load on the computer when simulating complex systems with high resolution. The high resolution achievable in finite-element programs cannot be achieved. Yet, unlike the finite-element programs, this method allows for the assessment of the mutual influences of thermal, electrical, and thermoelectric processes during design. If the accuracy of the results is unsatisfactory, more accurate calculation methods should be used instead.
Currently, the proposed method is in the initial stage of its development and requires further research. In particular, a more advanced method of convective heat transfer modeling needs to be developed, as well as a mesh optimization technique.