Peridynamic Formulation for Coupled Thermoelectric Phenomena

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.


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][3][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][6][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][11][12][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][15][16][17][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.

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 X A , where the capital superscript A denotes the material point index, it interacts with all its neighboring material points within a distance  called horizon radius.   depicts the family of X A [20].X B denotes the neighboring particles; the vector  AB = x B − x A is the bond vector from material point A to B. 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 X A 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 X A and all of its neighboring material points X B .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  AB and is defined as where Θ is the temperature.The temperature change between X B and X A 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  = 0, we obtain If () 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 X A .The virtual variation in  is given by Substituting ( 4) and ( 5) into (3), we obtained By defining the heat flow state as  = −∇(), we write the nontransient form of the balance of energy and electric charge without the source term as ( (X A , ) ⟨ AB ⟩ −  (X B , ) ⟨ BA ⟩)   = 0, (7a) 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  0 and  0 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 ( → 0) [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 X A is given, respectively, by where Φ⟨⟩ is PD electrical potential scalar state in the bond , and it is given by Φ(X A , )⟨ AB ⟩ = Φ(X B , ) − Φ(X A , ) and  is the shape tensor for thermal field and  EL 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.

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 as where  is the internal energy per unit mass and it is given by  =  V  =  V Θ,  is the temperature,  V 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,  0 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 as In 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 as For thermoelectric phenomena the heat flux  and the charge flux  can be obtained from the generalized Fourier's and Ohm's laws, respectively, as where  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 as 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: (24) 2.3.Bond Base Peridynamic (PD) Formulation for Thermoelectric Phenomena.In a bond-based peridynamic model, material point X A can interact with all neighboring material points X B 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 X A .κ 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: (26)

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 X A and X B 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 1 if there is an interaction among X A and X B 0 if there is no interaction among X A and X B . (28)

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, as where , ℎ, and  are the cross-sectional area and thickness and horizon radius, respectively.

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, A signifies our point of interest, and B signifies the points within the horizon of A.   is subdomain volume associated with the material point X B .For the purpose of PD simulation, the forward difference scheme is used.When forward differencing is employed, the following equation is solved:

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  = 0)).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.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 ( = 3.015Δ).

Initial
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 ,  V() , and  () are similar to that of example 1.The remaining material properties for bismuth telluride (Bi 2 Te 3 ) 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.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.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) 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 10 −4 sec.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.
The temperature inside the plate was initially 0 ∘ C and once 100 ∘ C 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 /2.
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 P and N 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 (  = −  ,   > 0) [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.,  = 0.05 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].
The domain is discretized in to 44 material points in the  and 20 in the  with time step of 10 −4 sec.
Figure 10 shows the maximum and minimum temperatures inside the thermopair along  = 0.7 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 (P and N).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  = 0)).To demonstrate effectiveness of PD theory in handling discontinuities, we considered a 2D thermoelectric plate of length  = 20 mm and width  = 20 mm with an insulated crack of length 2 = 10 mm, as shown in Figure 11.The plate is subjected to a temperature of −100 ∘ C at the bottom (cold junction) and 100 ∘ C 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

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 −35 ∘ C 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.

Figure 1 :
Figure 1: Broken thermal bond among material points X A and X B when intersecting insulated crack.

LFigure 2 :
Figure 2: Domain discretization of the real region.

Table 1 :
Geometric parameters and material properties.

Table 2 :
Material properties and geometric parameters.

Table 3 :
Geometric parameters and material properties.