Abstract
Modeling of heat and electrical current flow simultaneously in thermoelectric convertor using classical theories do not consider the influence of defects in the material. This is because traditional methods are developed based on partial differential equations (PDEs) and lead to infinite fluxes at the discontinuities. The usual way of solving such PDEs is by using numerical technique, like Finite Element Method (FEM). Although FEM is robust and versatile, it is not suitable to model evolving discontinuities. To avoid such shortcomings, we propose the concept of peridynamic theory to derive the balance of energy and charge equations in the coupled thermoelectric phenomena. Therefore, this paper presents the transport of heat and charge in thermoelectric material in the framework of peridynamic (PD) theory. To illustrate the reliability of the PD formulation, numerical examples are presented and results are compared with those from literature, analytical solutions, or finite element solutions.
1. Introduction
Depletion of fossil fuel, the global need for sustainable energy, and climate change due to the overutilization of fossil fuels are the major burning issues that speed up the research of new energy materials that are sustainable and environment friendly. From this point of view, thermoelectric materials become promising candidates due to their compactness, simplicity, scalability, reliability, environmental friendliness, portability, and excellent noiseless and maintenance-free properties. Furthermore, thermoelectric materials (TEMs) convert temperature differences into electric voltage; hence, it can theoretically convert waste heat to electricity or electricity to heat. We call the first phenomenon the Seebeck effect, where electrical voltage is developed when there is a gradient of temperature between the terminals and was discovered by Seebeck in 1821 while the second is just converse effect of the first one and it was discovered by Peltier in 1834 [1].
Despite the above advantages, TEMs are not without challenges due to their low energy conversion efficiency. Hence, substantial progress is underway to improve the performance of TEMs [2–4]. The efficiency of converting heat into electricity depends on temperature difference over which the device is operating, the geometrical design of thermoelectric device, and the figure of merit, which in turn depends on electrical conductivity, Seebeck coefficient, and thermal conductivity [5–7].
In general, TEMs are expected to operate in a very large temperature difference. As a result of this the device is subjected to huge amount of thermal stresses and strains for long periods of time resulting cracks at the interface. Hence, TEMs must be designed to withstand a large number of thermal cycles that cause mechanical fatigue [7]. Therefore, the life cycle estimation of thermoelectric materials due to fatigue crack is quite critical from the design point of view [8]. Finite element method has been widely used in the past for the life prediction of materials based on the classical theories. But these theories disregard the influence of defects in the material due to their divergence based governing equations. To overcome the shortcomings encountered by classical theories, Silling [9] introduced a new theory called peridynamics (PD) theory. Unlike classical theories that require some special consideration to define spatial derivatives of field variables at the discontinuities [10–13], peridynamics theory does not require any additional damage theory to capture discontinuities. Therefore, in this study the coupled transport of heat and charge in thermoelectric device will be presented using peridynamic theory. This is because the classical models [14–18] lead to infinite heat and electrical fluxes at the discontinuities. The theoretical framework developed here leads to a transient heat and charge transport models that do not contain spatial derivative and therefore are suitable to easily treat problems with discontinuities in the electric and temperature fields.
2. Materials and Methods
2.1. Review on State-Based Peridynamic Theory for Thermal Field
The term PD was originally proposed by [9] at Sandia National Laboratories in the late nineties to deal with the emergence and propagation of discontinuity in solids. According to Silling, the pairwise interaction is restricted to elastic materials with linear isotropic property and the effective Poisson ratio of 1/4. To avoid the limitations above, Silling et al. [19] proposed a generalized state-based PD approach in which the force in a bond is dependent on the cumulative deformation of all participating bonds connected with the material points within the horizon.
In this paper we follow the classical notion of continuum mechanics in Cartesian coordinates. Components that are defined in the undeformed configuration are represented by upper case indices and those quantities defined in the deformed configuration by lower case variables with subscripts .
For a material point , where the capital superscript denotes the material point index, it interacts with all its neighboring material points within a distance called horizon radius. depicts the family of [20]. denotes the neighboring particles; the vector is the bond vector from material point to . In the case of thermal problems, the temperature state at a material point is not measured by the local temperatures but by the temperature scalar state function .
Based on state-based PD theory [20, 21], the thermal potential of a material point is given by the functional denoted by and it is a function of its temperature state, which in turn depends on the temperature differences of and all of its neighboring material points . Therefore, the total thermal potential at a certain material point can be defined as a nonlocal integration over all bond vectors within the horizon , where is the thermal potential function and is the temperature scalar state operated on and is defined aswhere is the temperature.
The temperature change between and at any time is described by the temperature scalar state . Generally depends on temperature change associated with each interaction of a particular material point.
Invoking the stationary value of the thermal potential functional by setting , we obtainIf is a virtual variation in then is given by where represents the Fréchet differentiation of the thermal potential and is volume associated with the material point . The virtual variation in is given bySubstituting (4) and (5) into (3), we obtainedBy defining the heat flow state as , we write the nontransient form of the balance of energy and electric charge without the source term as where is current flow state.
Equations (7a) and (7b) may be considered as the classical complements of the steady state balance of energy and the balance of electric charge equations ((8a) and (8b)) without heat and charge sources, respectively,where and are the heat and charge fluxes in reference configuration, respectively. It has been proven that the classical equations ((8a) and (8b)) and the peridynamic equations ((7a) and (7b)) are the same in the limit as the horizon radius approaches zero [22, 23]. For peridynamic heat and electrical conduction phenomena, the discrete form of the nonlocal temperature and electrical potential gradients for a certain material point is given, respectively, bywhere is PD electrical potential scalar state in the bond , and it is given by and is the shape tensor for thermal field and is the shape tensor of electric field. Equations (9a) and (9b) are crucial steps in developing the nonlocal temperature and electrical potential states at a certain material point.
2.2. Peridynamic Formulation for Thermoelectric Materials
Equations (7a) and (7b) are the balance of energy and charge equations for steady state problems without considering the heat and charge source terms. The energy and charge balance equations including the heat and charge source term in the reference configuration can be described aswhere is the internal energy per unit mass and it is given by , is the temperature, is the specific heat capacity at a constant volume, is the Cauchy heat flux in current configuration, is the heat source term per unit mass per unit time, is the mass density, is the electrical charge of the carrier, is the electrical charge flux, is the charge source term per unit mass per unit time, is deformation gradient, and is the Jacobian.
By relating the LHS of (7a) and (7b) and (8a) and (8b) and replacing the nonlocal integral with finite sum the PD governing equations may be expressed asIn thermoelectric phenomena, the coupling between the gradients of temperature and electrical potential gives rise to the thermoelectric effect of Seebeck, Peltier and Thomson. For the steady state condition the energy and charge balance equations for thermoelectric phenomena are given, respectively, by For transient condition, (14a) and (14b) are modified asFor thermoelectric phenomena the heat flux and the charge flux can be obtained from the generalized Fourier’s and Ohm’s laws, respectively, aswhere is thermal conductivity of the material, is the electric conductivity of the material, is electric potential, is temperature, and is Seebeck coefficient. Therefore, (16) may be expressed as follows:For transient conditions including the heat and charge source term the energy and charge balance equations for thermoelectric phenomena may be written in the current configuration asKnowing To carry out the peridynamics formulation of thermoelectricity, we first rewrite (19a) and (19b), in the reference configuration as follows:Therefore, the PD governing equation for thermoelectric phenomena may be expressed as follows: where the first, the second, and the third summations are equivalent to the following classical equations in the same order:Equation (21) is the general form of the state-based PD heat transfer equation for thermoelectric phenomena.
Similarly, the state-based PD equation for the transport of charge in thermoelectric phenomena may be obtained as follows: where the first and the second summations are equivalent to the following classical equations in the same order:
2.3. Bond Base Peridynamic (PD) Formulation for Thermoelectric Phenomena
In a bond-based peridynamic model, material point can interact with all neighboring material points in its horizon in a pairwise manner. The change in temperature at the two end points of a bond is assumed to cause the heat to flow along the axis of the bond only. When material points interact in a pairwise manner and are restricted to a specified neighborhood through a bond, (21) may be reduced as follows: and are microconductivity of the electrical and thermal bonds, respectively, and is the horizon volume of material point centered at .
and are the PD conductivity of thermal and electrical bonds, respectively, between material points and . When we decouple (25) by disregarding the Seebeck effect, we obtain the PD heat conduction equation of [20] as follows:
2.3.1. Damage Modeling
A most common damage in engineering is the presence of cracks in the material. In the case of heat conduction phenomena, the presence of cracks in a body will terminate the heat conduction process completely or partially. Cracks that terminate heat transfer are said to be insulated cracks. In our peridynamic model, the crack surfaces are created by terminating the microthermal potentials between material points and when crossing the insulated crack C-D as shown in Figure 1. This can be done by removing the peridynamic heat flow densities of these particular sets of material points from the PD heat conduction equation (26). Hence, (26) may be written by incorporating the history dependent scalar value function [23] as follows:where
2.3.2. Linking Peridynamic Properties with That of the Classical Counterparts
In order to create a connection between the PD properties and the classical material properties we directly borrow expressions from [20]. Thus, for one- and two-dimensional analysis, the PD thermal and electrical microconductivities are expressed, respectively, aswhere , , and are the cross-sectional area and thickness and horizon radius, respectively.
3. Results and Discussions
3.1. Numerical Procedures
To solve the PD heat conduction equation for thermoelectric phenomena, numerical techniques are used. The domain is discretized using equally spaced grids for simplicity. Equation (25) can be numerically solved by replacing the nonlocal integral equation with finite sum as follows:where is the number time step, signifies our point of interest, and signifies the points within the horizon of .
is subdomain volume associated with the material point . For the purpose of PD simulation, the forward difference scheme is used. When forward differencing is employed, the following equation is solved:
3.2. Validation of PD Theory
To show the applicability of the proposed peridynamics approach, we implemented the bond-based peridynamics approach. Six numerical examples are presented; in the first example the uncoupled Seebeck effect for one-dimensional bar has been carried out and the results from PD simulation are compared with the analytical solution. In the second example the coupled Seebeck effect for the coupled thermoelectric phenomena has been carried out and the results from PD simulation are compared with results from the literature. The third example demonstrates the effect of Joule heating in thermoelectric materials for the case of nonsymmetric temperature boundaries. The fourth example illustrates two-dimensional uncoupled Seebeck effect for symmetric temperature boundaries. The fifth example demonstrates the coupled Peltier effect for thermoelectric phenomena. Finally, we presented the capability of PD theory in handling discontinuities.
Example 1 (heat conduction of one-dimensional thermoelectric bar (uncoupled case when )). Consider a discretized one-dimensional thermoelectric bar of length as shown in Figure 2 with material properties and geometric parameters listed in Table 1. The domain is discretized into 100 material points and the spacing (Δ) between them is 0.0001 m with time step of 0.01 sec.
Initial Conditions. , .
Boundary Conditions. , .
To apply the boundary condition through the fictitious regions, we need to add extra material points for the boundaries [20]. As shown in Figure 3, three material points are added on the left and on the right (pink cells with dashed lines in Figure 3). As recommended by [20], the size of the fictitious region may be chosen as the size of the horizon (δ), which is equal to 3.015 times the spacing ().
in Figure 3 is the fictitious boundary [20]. With reference to Figure 4, temperature variation results from the PD transient heat conduction and analytical steady state solution are in good agreement. It is also noticed from Figure 4 that temperature increases towards the right boundary as expected.
Example 2 (heat conduction of one-dimensional thermoelectric bar (coupled Seebeck effect)). In this example, the domain is discretized into 20 material points and spacing between them is 0.0005 m with time step of 10 μsec. The boundary conditions , , and are similar to that of example 1. The remaining material properties for bismuth telluride (Bi2Te3) are listed in Table 2 [18].
As shown in Figure 5, the PD transient heat conduction results for temperature dependent material properties are compared with results obtained from [18] with constant material properties and found to be in good agreement.
Example 3 (heat conduction of one-dimensional thermoelectric bar with constant electric flux). In this example, we considered two boundary cases (symmetric and nonsymmetric boundaries). The material properties and the geometrical parameters are similar to those in Example 2. The constant electric flux is induced by a voltage supply () at left boundary. The domain is discretized into 50 material points and spacing between them is 0.03048 mm with time step of 0.0001 sec.
Boundary Conditions. , , , .
Figure 6 shows the PD and steady state analytical result of thermoelectric bar subjected to nonsymmetric boundary. We also observed from Figure 6 that the temperature value increases at the center due to Joule heating effect. This effect is one of the irreversible thermoelectric effects and its magnitude depends basically on the current and electrical conductivity of the material.
Example 4 (heat conduction of two-dimensional plate). To demonstrate capability of PD theory we extend the uncoupled one-dimensional problem to two-dimensional problems. The material properties and geometrical parameters are indicated in Table 3. In this example again we considered two boundary cases, symmetric and nonsymmetric boundaries. In case of nonsymmetric boundary, the following boundary and initial conditions are considered.
Case 1 (nonsymmetric boundary)
Initial Conditions. , , .
Boundary Conditions. , .
The domain here is discretized into 20 material points in the and 20 in the . Spacing between them is 0.1 cm with time step of . Figure 7 shows two-dimensional temperature variations from PD and analytical solutions. It is observed that results from the two solutions are in good agreement and also the temperature rises towards the right boundary as we expect.
Case 2 (symmetric boundary). When we consider symmetric boundary, the following boundary condition is used.
Boundary Conditions. , .
The temperature inside the plate was initially and once temperature is imposed on the boundaries, the temperature inside the plate increases starting from the ends. As shown in Figure 8, the temperature increases with time and attains its maximum value at .
Example 5 (thermoelectric cooler). When an external voltage is supplied at one end of the thermopair, one of the junctions liberates heat and the other absorbs heat. Therefore, by supplying an electric power we can get a cooling or heating effect by using thermoelectric cooler. In this example, the and type thermoelements are joined to form a thermopair as shown in Figure 9. The thermal and electric conductivities of the thermoelements are considered to be the same but their Seebeck coefficient α is different [16]. Our main objective in this example is to find the temperature distribution at the junction. We technically consider the problem as one-dimensional case as far as this example is concerned.
We considered three different cases based on the voltage supply at the left boundary (i.e., v, 0.15 v, and 0.25 v). The voltage at the right boundary is assumed to be 0 v for all the three cases. Similarly, the temperature boundaries on all sides of the domain are assumed to be 7°C. The material properties and geometrical parameters are listed in Table 4 [16].
Boundary Conditions. , , , .
The domain is discretized in to 44 material points in the and 20 in the with time step of .
Figure 10 shows the maximum and minimum temperatures inside the thermopair along mm for the three voltage sources (0.05 v, 0.15 v, and 0.25 v). We observed that the influence of Joule effect increases as we move from left to right and it reaches its maximum value at the middle of each thermoelement ( and ). We also noticed that there is a considerable temperature drop at the junction due to Peltier effect. The result pattern in this study follows a similar pattern with those in [16], but the numerical values slightly differ. This is due to the fact that we reduced the 3D problem in [16] to 2D problem.
Example 6 (heat conduction of two-dimensional thermoelectric plate with insulated crack (when )). To demonstrate effectiveness of PD theory in handling discontinuities, we considered a 2D thermoelectric plate of length mm and width mm with an insulated crack of length , as shown in Figure 11. The plate is subjected to a temperature of at the bottom (cold junction) and at the top (hot junction), by keeping the right and left boundaries insulated. The material properties are taken from Table 1.
The domain here is discretized into 100 material points in the and 100 in the . Spacing between them is 0.02 cm with time step of 0.01 sec. The same problem has been solved using ABAQUS to show the validity of our PD result for verification purpose. In this finite element model the domain is discretized into 9900 linear quadrilateral elements. Figures 12(a) and 12(b) show the temperature distributions from PD and ABAQUS solutions, respectively. It is observed from Figures 12(a) and 12(b) that the temperature contour lines near the insulated crack deviated from straight horizontal line as we get closer to the crack. This happened due to the fact that the flow of heat is affected by the presence of the crack. It is also observed that results from the two solutions are in a very good agreement.
(a)
(b)
4. Conclusion
In this paper we presented the peridynamic formulation of coupled heat and charge transport in thermoelectric media. The local gradients of temperature and electrical potentials are replaced by the functional integral of the temperature and electrical potential fields that are valid even in case of evolving discontinuities. In this study the material is assumed to be isotropic and homogeneous. Based on these assumptions, we derived the peridynamic balance and constitutive equations in thermoelectric medium. Furthermore, we developed the peridynamic equations for the conservation of energy and charge for coupled thermoelectric phenomena. To demonstrate the capability of the proposed PD formulation we solved six examples and compared the results with analytical solutions and results from the literature or from finite element solutions. In the first example we considered 1D heat conduction for the uncoupled thermoelectric bar. In this example, we compared our result with analytical solution and found that they are in close agreement. The second example solved a 1D heat conduction problem for the couple thermoelectric phenomena. It was observed in this example that our result and result from the literature were in close agreement. The third and fourth examples presented 1D heat conduction with constant electric flux and heat conduction in 2D plate, respectively. Results from PD solution were compared with the analytical solution and they are found to be in good agreement. Regarding example five, interesting results are obtained at the junction of thermopair as we expect. We observed that temperature attained its maximum value at the center of each thermoelement (P and N) due to Joule heating effect. We also noticed that there is a considerable temperature drop of at the junction due to Peltier effect. In the last example, we demonstrated the ability of PD in handling discontinuities. In this example, we modeled 2D thermoelectric plate with insulated crack. Results from PD and ABAQUS were compared and found to be in good agreement. Since the derived PD equations are valid whether there are discontinuities or not in the domain of interest, PD has a huge potential in solving problems related to thermal fatigue. In the near future we will employ peridynamics to include elastic field in addition to thermal and electric fields.
Disclosure
Part of this work has been poster presented during the 35th International Conference & the 1st, Asian Conference on Thermoelectrics (ICT/ACT 2016) in Wuhan.
Competing Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
This work is financially supported by the National Program on Key Basic Research Project (973 Program, no. 2013CB632505), National Natural Science Foundation of China (no. 10672127), National Natural Science Foundation of China (U12301013), and Chinese Scholarship Council (CSC).