Non-linear finite element formulation applied to thermoelectric materials under hyperbolic heat conduction model

Comput. Methods Appl. Mech. Engrg. 213–216 (2012) 93–103

Contents lists available at SciVerse ScienceDirect

Comput. Methods Appl. Mech. Engrg. journal homepage: www.elsevier.com/locate/cma

Non-linear finite element formulation applied to thermoelectric materials under hyperbolic heat conduction model R. Palma a,⇑, J.L. Pérez-Aparicio a, R.L. Taylor b a b

Mecánica de Medios Continuos y Teoría de Estructuras, Universidad Politécnica de Valencia, Spain Department of Civil and Environmental Engineering, University of California at Berkeley, CA, USA

a r t i c l e

i n f o

Article history: Received 9 March 2011 Received in revised form 12 August 2011 Accepted 14 November 2011 Available online 25 November 2011 Keywords: Thermoelectrics Dynamics Non-linear finite element Second sound HHT algorithm Newmark-b algorithm

a b s t r a c t In the present work, a three-dimensional, dynamic and non-linear finite element to simulate thermoelectric behavior under a hyperbolic heat conduction model is presented. The transport equations, which couple electric and thermal energies by the Seebeck, Peltier and Thomson effects, are analytically obtained through extended non-equilibrium thermodynamics, since the local equilibrium hypothesis is not valid under the hyperbolic model. In addition, unidimensional analytical solutions are obtained to validate the finite element formulation. Numerically, isoparametric eight-node elements with two degrees of freedom (voltage and temperature) per node are used. Non-linearities due to the temperature-dependence on the transport properties and the Joule effects are addressed with the Newton–Raphson algorithm. For the dynamic problem, HHT and Newmark-b algorithms are compared to obtain accurate results, since numerical oscillations (Gibbs phenomena) are present when the initial boundary conditions are discontinuous. The last algorithm, which is regularized by relating time steps and element sizes, provides the best results. Finally, the finite element implementation is validated, comparing the analytical and the numerical solutions, and a three-dimensional example is presented. Ó 2011 Elsevier B.V. All rights reserved.

1. Introduction Thermoelectric materials couple electric and thermal energies by means of three separated transport effects: Seebeck, Peltier and Thomson. In addition, the Ohm and Fourier laws that are inherent to electric and thermal fluxes are also present. Thermoelectric materials are used as heat pumps (heating and cooling) and generators, see [1], and in the last decade miniaturized thermoelectric devices and high-frequency processes are increasingly applied to cool micro-electronic devices, see [2]. The classical Fourier law leads to a parabolic heat propagation problem incorrect from a physical point of view: especially applied to micro-devices and under rapid transient effects such as micropulses. Cattaneo [3] proposed a modification to the Fourier law that leads to a hyperbolic heat propagation problem, avoiding the classical paradox related with the infinite temperature velocity of propagation. The modification, named sometimes ‘‘second sound’’, is made through the introduction of a relaxation time s, an empirical parameter defined as the time-interval between two successive collisions at the microscopic level of either holes or electrons. From a theoretical point of view, the classical non-equilibrium thermodynamics [4] cannot lead to hyperbolic propagation since the local equilibrium hypothesis is not valid. For this reason, ex⇑ Corresponding author. E-mail address: [email protected] (R. Palma). 0045-7825/$ - see front matter Ó 2011 Elsevier B.V. All rights reserved. doi:10.1016/j.cma.2011.11.011

tended non-equilibrium thermodynamics [4,5] must be applied to obtain the correct transport equations. Several finite element (FE) formulations applied to the thermoelectric problem using the Fourier (parabolic) model, see [6,7] have been published. Analytically, in [8] the Laplace transform technique was used to solve problems formulated using the Cattaneo (hyperbolic) model. For heat propagation (without thermoelectric coupling) analytical investigations were performed using the Laplace transform technique under the Cattaneo model [9–11]. In [12] a mixed FE formulation using the Crank–Nicolson scheme for timeintegration was presented and in [13] a finite difference scheme was developed to compare implicit and explicit Euler algorithm solutions. Several FE works in which the Cattaneo model was used also have been published. For example, in [14,15] FE’s were formulated to simulate convection–diffusion (Fick’s law) and thermomechanical coupling, respectively. The aim of the current work is to develop a three-dimensional (3D), transient, non-linear FE formulation to solve general thermoelectric problems with the hyperbolic model. For this purpose, a non-equilibrium entropy depending on temperature and heat flux is defined and the transport equations are formulated by applying the Onsager relations. Transport and balance equations are then transformed into a matrix form and implemented into the research FE code FEAP [16]. The non-linearities due to: (i) Joule effect, (ii) temperature-dependence of transport properties and (iii) thermoelectric coupling between electric and heat fluxes, are solved with

94

R. Palma et al. / Comput. Methods Appl. Mech. Engrg. 213–216 (2012) 93–103

the Newton–Raphson (N–R) algorithm. Two discrete algorithms for transient solutions, namely, HHT or a [17] and Newmark-b (N-b) are compared. Both are regularized by means of a linear relationship between time steps and element lengths. In order to validate the FE behavior, four unidimensional (1D) analytical solutions are compared with the numerical results. When the initial boundary conditions are discontinuous, strong oscillations (Gibbs phenom and c  ena) are observed and the N-b algorithm, with modified b parameters, provides results that are consistent and accurate with respect to the analytical solutions. Finally, a 3D case is simulated, obtaining several conclusions for the design and optimization of thermoelectric devices.

Neglecting qX, assumption valid for most thermoelectric applications [19], the charge and energy balances are

The extended non-equilibrium thermodynamics theory assumes the existence of a non-equilibrium entropy s that depends on state variables and on fluxes. Since the electric energy propagates much faster than the thermal one, the hyperbolic model is applied only for the thermal field, and the entropy balance is given by Jou and Lebon [4] adding the Joule irreversibility

qs_ þ r 

2.1. Balance laws Consider a thermodynamic universe composed of a system of domain X, boundary C and its surroundings, see Fig. 1. Assuming the validity of a continuum hypothesis, the balance laws of continuum physics are those of charge and energy. The first is obtained from a Maxwell equation, see [18]

q_ X ¼ r  j;

ð1Þ

where a dot () denotes a time derivative, j is electric flux and qX is charge density. Energy balance is obtained by applying the first law of thermodynamics

qu_ ¼ r  1 þ j  E;

ð3Þ

2.2. Transport equations

2. Thermodynamic formulation The aim of this section is develop balance and transport equations for the hyperbolic thermoelectric problem. As mentioned above, the hyperbolic problem cannot be obtained from classical non-equilibrium thermodynamics. This is due to the lack of validity of the local equilibrium hypothesis that states the equality of the local relations between thermodynamic quantities in a system in and out of equilibrium [5]. This validity is closely related with the Deborah number De = s/sM < 1, where sM is the macroscopic experimental time (effect duration) and s the time interval between two successive collisions between particles. Large Deborah numbers occur if: s is large, for materials such as polymers, and/or sM is short, for very fast phenomena such as ultrasound propagation. The extended non-equilibrium thermodynamics [5] must be used for transport equations when De P 1.

qu_ ¼ r  1 þ j  E:

r  j ¼ 0;

1 T

¼

    j  rV 1 s  2 1_ ; þ1 r T T jT

ð4Þ

where T, V, j denote temperature, voltage and thermal conductivity, respectively. From the general form of the entropy balance, the entropy flux js and the rate of entropy production rs are

js ¼

1 T

;

rs ¼

    j  rV 1 s þ1 r  2 1_ : T T jT

ð5Þ

The classical form of the entropy production for the thermoelectric problem, see [20], is recovered if s = 0: the hyperbolic model adds a term to the entropy production. Finally, the transport equations are obtained by expressing the entropy production in the form from [5] and by satisfying the Onsager–Casimir reciprocal relations to give

j ¼ crV  acrT; 1 þ s1_ ¼ jrT þ aTj;

ð6Þ

where a, c denote the Seebeck coefficient and thermal conductivity, respectively. These phenomenological coefficients are denominated thermoelectric properties in the following. Again if s = 0 in (6)-bottom, the classical thermoelectric transport equations are recovered. Although linear relationships between fluxes and gradients have been used, the thermoelectric properties depend strongly on T (material non-linearity), accordingly in the present work

aðTÞ ¼ 1:98  104 þ 3:35  107 T  7:52  1010 T 2 ; cðTÞ ¼ 1:09  105  5:59  102 T þ 2:49T 2 ; jðTÞ ¼ 1:66  3:58  103 T þ 3:19  105 T 2

ð7Þ

in which dimensions are: a [V/K], c [A/mV], j [W/mK] but T is introduced in Celsius (°C) degrees, see [21].

ð2Þ

_ 1 and E denote mass density, rate of internal energy where q, u, density, heat flux and electric field, respectively. The last term on the right-hand represents the rate of electromagnetic energy converted into thermal by the irreversible Joule effect.

2.3. Strong form From the previous equations and taking into account the coupled nature of the thermoelectric problem, the strong form is formed as a system of two coupled partial differential equations. The electric strong form is obtained by replacing (6)-top into (3)left and considering the T-dependence of the thermoelectric properties (7)

cðTÞr2 V þ @ T fcðTÞgrT  rV ¼ ½cðTÞ@ T faðTÞg þ aðTÞ@ T fcðTÞg þ aðTÞcðTÞr2 T: ð8Þ Here and in the following, @ T implies partial derivation with respect to the subscript or subscripts, in this equation the temperature T. To obtain analytical solutions, from (6)-top we deduct

rV ¼ aðTÞrT  Fig. 1. Energy balance in a domain X: the three effects, reversible couplings (arrows without border) and irreversible thermal and electric fluxes (arrows with border).

j

cðTÞ

;

ð9Þ

where the first term in the right-hand represents the thermal energy converted into electric by the Seebeck effect. Due to restrictions

R. Palma et al. / Comput. Methods Appl. Mech. Engrg. 213–216 (2012) 93–103

in analytical solutions (not in the numerical ones), in this article we present cases in which only a is variable with T. The values of j, c will be evaluated at a certain average T from (7). The thermal strong form is obtained substituting (6)-bottom into (3)-right and using (3)-left. Furthermore, the relations u = cT (c, specific capacity), E = rV are introduced

j



95

R

rdV  jdX þ C dVjc dC ¼ 0; R R X rdT  qdX þ X dTj  rVdX þ X dT s@ t fj  rV gdX

RX þ þ

R

RX C

€ Xþ dT sqcTd

R

X

_ X dT qcTd

ð13Þ

dTqc dC ¼ 0:

3.2. Discretization

2

sqcT€ þ qcT_ ¼ jr2 T þ  th rT  j  s@ t fj  rVg; c



R

ð10Þ The continuum domain X is discretized into elements Xe, see nel Fig. 2. Therefore, X ¼ [e¼i Xe where nel is the total number of elements, each one delimited by npe = 8 nodes. In turn, each node relates to two degrees of freedom (voltage, temperature) per node. For an element e, degrees of freedom and spatial coordinates x are approximated using 3D isoparametric shape functions Ne

where th = T[@ a(T)/@T] is the Thomson coefficient. This expression is an hyperbolic equation, notice the first term on the left-hand. The right-hand terms represent several effects: the first classical heat conduction (Fourier), the second Joule, the third Thomson and the last dissipation related to the hyperbolic model. Fig. 1 schematizes the energy balance of the thermodynamics formulation. Electric and thermal energies are represented left and right, respectively; sources top and the flux bottom. The electric energy is irreversibly converted into thermal by the Joule effect, therefore the related arrow points towards the thermal energy side. Inside the domain, electric energy is reversibly converted into thermal and vice versa by three separated effects: Seebeck (thermal into electric), Peltier (electric into thermal) and Thomson in both ways.

n o n o where Ve ¼ V 1e ; . . . ; V ne pe , Te ¼ T 1e ; . . . ; T ne pe . Be is the discretized gradient matrix of the element composed of npe submatrices with the form:

3. FE formulation

BtA ¼

The aim of this section is to develop a variational formulation within the FE framework to permit the implementation of the thermoelectric problem into a computer code. For this purpose, consider a thermoelectric body of domain X and boundary C governed by balance and transport equations obtained in the previous section. To facilitate the FE implementation, again the hyperbolic flux s1_ from the thermal transport Eq. (6)-bottom is introduced in the energy balance (3)-right, with the result that the expression of the pure parabolic heat flux q is recovered in the transport equation

sqcT€ þ qcT_ ¼ r  q  j  rV  s@ t fj  rVg q ¼ jrT þ aTj

ð11Þ

V  V h ¼ Ne V e ;

T  T h ¼ Ne Te ;

h

rV  rV ¼ Be Ve ; T_  T_ h ¼ Ne T_ e ;

rT  rT h ¼ Be Te ;

€e ; T€  T€ h ¼ Ne T

ð14Þ

x  xh ¼ Ne xe ;

h

@N A @x1

@N A @x2

@N A @x3

i ;

ð15Þ

where ()t denotes matrix transpose. Similarly, the variations are discretized as

dV  dV h ¼ Ne dVe ; dT  dT h ¼ Ne dTe ; rdV  rdV h ¼ Be dVe ; rdT  rdT h ¼ Be dTe :

ð16Þ

Introducing these discretizations into the transport equations (6)top and (11)-bottom

je ¼ cBe Ve  acBe Te ; qe ¼ jBe Te  acNe Te Be Ve  a2 cNe Be T2e :

ð17Þ

Introducing now (14)–(16) into (13) the Galerkin forms are obtained

The governing equations are complemented by Dirichlet (on boundaries CV, CT) and Neumann (on Cj, Cq) boundary conditions for electric and thermal fields

V¼V

on CV ;

j  n ¼ jc

on Cj ;

T¼T

on CT ;

q  n ¼ qc

on Cq ;

where V, T, jc and qc are prescribed voltage, temperature, electric flux and thermal flux, respectively, and the outward normal is denoted by n. Boundaries satisfy the conditions CV [ Cj = C, CV \ Cj = 0 and CT [ Cq = C, CT \ Cq = 0. 3.1. Weak forms According to standard variational methods [22], the weak forms are obtained by multiplying the balance Eq. (3)-left and (11)-top by variations of the degrees of freedom, dV, dT and integrating over X

R R

dV r  jdX ¼ 0; R R _ X dT r  qdX þ X dTj  rVdX þ X dT qc TdX R R € þ X dT s@ t fj  rVg dX þ X dT sqcTdX ¼ 0: X

ð12Þ

The s and q-terms are new with respect to the formulation developed in [7]. The divergence theorem is now applied to the first term of both equations, and Dirichlet boundary conditions are introduced

Fig. 2. Domain discretization in finite elements Xe. Isoparametric three-dimensional and interface flux elements. Two degrees of freedom per node.

96

R. Palma et al. / Comput. Methods Appl. Mech. Engrg. 213–216 (2012) 93–103

 

R R

R

Xe

Bte je dXe þ

Xe

Bte qe dXe þ

þ þ

R R

Nte

Ce

R

Nte jc dCe ¼ 0;

V K VV AB ¼ @ V B fRA g ¼

R

Nte je Be Ve dXe þ Xe sNte @ t fje Be Ve gdXe R qcNe T€ e dXe þ Xe Nte qcNe T_ e dXe Xe

Xe

s

Ce

Nte qc dCe ¼ 0:

ð18Þ

3.3. Non-linear transient solution form In the present work a non-linear, transient hyperbolic problem is studied. As mentioned several non-linearities are present, and the transient behavior is relevant only for the thermal field, since qX in (1) is assumed to be zero. In order to solve the transient problem, the time interval of interest is divided into small increments, Dt, and the N-b or HHT method is introduced to replace time derivatives by discrete forms. Then for each time point the N–R method is used to solve the resulting non-linear algebraic problem. These transient algorithms have been widely used in structural and solid mechanics (see [22]) and will not be described in detail here. Basically, the assembled non-linear FE equations are written in a residual form R and linearized by

k @RA  RkA ¼  dgkA g 

ð19Þ

B

where A, B are the global numbering of two nodes, k the N–R iteration counter and gA derivatives of the degrees of freedom at node A: _ A ¼ fV_ A ; T_ A gt , and second value as UA = {VA, TA}t, first derivative as U t € € derivative as UA ¼ f0; T A g . The algorithms for time integration are taken as

k @RA   ¼ c1 KAB þ c2 CAB þ c3 MAB : @gB 

ð20Þ

 c  for In Table 1 the parameters c1, c2 and c3 are given in terms of b,  for the HHT algorithm. the N-b algorithm and a The consistent tangent, capacity and mass matrices, are derived for each iteration

KAB ¼ 

@RA ; @UB

CAB ¼ 

@RA ; @ U_ B

MAB ¼ 

@RA : €B @U

ð21Þ

Finally, the solution is updated by gkþ1 ¼ gkB þ dgkB . Note that to enB sure a correct derivation of the tangent matrix, the N–R should exhibit a quadratic asymptotic rate of convergence. From (18), the assembled residuals for each field at each global node are formulated as

RVA ¼  RTA

R

BtA jA dX þ

R

C N A jc dC; R t B q d X þ X A A X N A j A BA V A dX R R þ X NA s@ t fjA BA V A gdX þ X NA sqcT€ A dX R R þ NA qcT_ A dX þ NA q dC:

¼

R

X

X

C

V K VT AB ¼ @ T B fRA g ¼

K TV AB

c

From (21)-left and (22), the tangent stiffness matrices are

BtA @ V B fjgdX;

X

BtA @ T B fjgdX;

R

ð23Þ

Similar tangent matrices were obtained in [7] for the steady-state problem, situation for which both theories (parabolic and hyperbolic) are equivalent. Using now (17) and the chain rule

@ V B fjg ¼ cBB ; @ T B fjg ¼ acBB  @ T fagNB cBB T B @ T fcgðBB V B þ aBB T B ÞNB ;

ð24Þ

@ V B fqg ¼ aNB T B @ V B fjg; @ T B fqg ¼ aNB T B @ T B fjg þ ajB NB  jBB þ@ T fagNB T B jB NB  @ T fkgBB T B NB : From (21)-middle and (22), the capacity matrices are

n o V C VV ¼ 0; AB ¼ @ V_ B RA C TT AB

n o V C VT ¼ 0; AB ¼ @ T_ B RA

R NA @ tV_ B fjgBB V B dX  s X NA jB BB dX; R R ¼ @ T_ B fRTA g ¼ s X NA @ tT_ B fjgBB V B dX  X NA qcN B dX;

T C TV AB ¼ @ V_ B fRA g ¼ s

R

X

ð25Þ where the derivatives are again obtained applying the chain rule to the discretized transport Eq. (17)

@ tV_ B fjg ¼ cBB ;

@ tT_ B fjg ¼ acBB :

ð26Þ

As mentioned the symbol @ tV_ B implies second partial derivative with respect to the time and the temporal derivative of V. From (21)right and (22), the mass matrices are obtained by V MVV € fRA g ¼ 0; AB ¼ @ V B

V M VT AB ¼ @ T€ B fRA g ¼ 0;

T MTV € fRA g ¼ 0; AB ¼ @ V B R T MTT AB ¼ @ T€ B fRA g ¼  X N A sqcN B dX:

ð27Þ

The final assembled matrix problem is given by

"

c1 KTT þ c2 CTT þ c3 MTT



T V

ð22Þ

X

R

¼ X BtA @ V B fqgdX R  X NA @ V B fjgBB V B dX  X N A jA BB dX R R  X NA s@ V B fjgBB V_ B dX  X NA s@ t fjgBB dX; R t T K TT AB ¼ @ T B fRA g ¼ X BA @ T B fqgdX R R  X NA @ T B fjgBB V B dX  X NA s@ T B fjgBB V_ B dX: ¼ R

@ V B fRTA g

R

kþ1

c1 KTV þ c2 CTV

c1 KVT c1 KVV k

k T dT ¼ þ : V dV

#

dT dV

(

k ¼

RT RV

)k ;

ð28Þ The present formulation is implemented as a user element in the non-linear finite element code FEAP, see [16]. Each element uses a standard, in the finite element method (FEM) sense, eight-node isoparametric element with two degrees of freedom per node A: voltage VA and temperature TA. 4. Validations

Table 1 Tangent matrix parameters for Newmark and HHT algorithms. Method

Parameters c1

c2

c3

N-b

1

HHT

a

c  Dt b a c  Dt b

1  Dt 2 b 1  Dt 2 b

The aim of this section is to validate the FE developed in Section 3 by means of four cases for which analytical solutions are found in the literature or modified and developed here. These solutions are 1D simplifications of the coupled partial differential equations (9) and (10) and are summarized in Table 2. Due to the difficulty in developing analytical solutions, the T-dependence of the electric and thermal conductivities is not included in the validations. However, for the steady-state problem, this dependence has been verified in [21] using experimental results.

R. Palma et al. / Comput. Methods Appl. Mech. Engrg. 213–216 (2012) 93–103 Table 2 Unidimensional validation cases, simplifications and considered effects: F - Fourier, S Seebeck, J - Joule, Th - Thomson, ss - second sound. Case

Simplifications

Effects

I II III IV

s = 0; j = 0; a, c, j = cte s = 0; j – 0; a, c, j = cte s = 0; j – 0; a(T); c, j = cte s – 0; j – 0; a(T); c, j = cte

F, F, F, F,

S S, J S, J, Th S, J, Th, ss

For all cases, a p-type thermoelement (TE) device manufactured by MELCOR [23] with properties described by (7) is modeled. The TE is a parallelepiped of dimensions 1.4  1.4  1.14 [mm]. For a p-type TE, the electric flux direction is co-linear with the heat flux. The boundary conditions are discontinuous: initially the temperature is T = 0 everywhere. For t > 0 the boundary temperatures are set to Tc = 30 on the cold face and to Th = 50 °C on the hot face. The boundary voltage is set to V = 0 [V] on the cold face. Fig. 3 shows the geometric dimensions and the boundary conditions. For cases II, III and IV, an electric flux in the x3 direction is prescribed using the special interface element developed by the authors in [21], see Fig. 2. The applied intensity 5.2 [A] is an average of the ones specified by the manufacturer and corresponds to an electric flux j3 = 2.65  106 [A/m2]. The constant parameters for a, c, j are obtained from (7) using an average temperature Tm = (Th + Tc)/2. For all cases, the electric strong form is directly given by (9), however, the thermal strong form for each case must be obtained by introducing the simplifications from Table 2 into (10)

j@ xx fTg ¼ qc@ t fTg; j@ xx fTg ¼ qc@ t fTg  j2 =c; 2 III : j@ xx fTg ¼ qc@ t fTg  j =c þ t h j@ x fTg; 2 IV : j@ xx fTg ¼ qc@ t fTg  j =c þ t h j@ x fTg þs@ t fj@ x fVgg þ sqc@ tt fTg:

I:

II :

ð29Þ

The fourth and fifth terms of the right side of case IV (second line) represent the second sound effects: irreversibility, (5)-right and hyperbolicity, respectively. 4.1. Cases I–III For these cases, s = 0 and the thermal strong forms (29) are diffusive: second order parabolic partial differential equations. For I, the situation is linear and homogeneous; for II linear and nonhomogeneous and for III non-linear and non-homogeneous. Analytical solutions for the thermal field in cases I and II are given in [24], while the electrical fields can be calculated from (9), see the Appendix. For III, the thermal and electrical solutions are given in [25]. Numerical solutions are obtained using a structured and coarse mesh of 11 elements in the x3 direction. Only one element is used

97

in the x1 and x2 directions since the problem is fundamentally 1D. A time step of Dt = 0.1 [s] and the standard N-b parameters  ¼ 0:25, c  ¼ 0:5 are used. b Fig. 4 shows the voltage (left) and temperature (right) distributions along the x3 direction for case I (top), and for coincident cases II and III (bottom). Solutions at t = 0.05, 0.3, 5 [s] are represented: the analytical ones with lines and the FEM results with circles. For case I, the T distributions are quadratic in nature during the initial transient response due to the boundary conditions. At near steady-state (e.g., t = 5 [s]) they become nearly linear, since this problem is of the Laplace-type: the lack of electric current implies the absence of a Joule effect. The V distributions are proportional to T, see (9). For case II, near quadratic distributions appear for the steadystate due to the Joule effect: this is a Poisson-type problem. Now, V distributions are not proportional to T, since the potential drop increases due to electric energy being converted into thermal energy. Results for case III are very similar, since the Thomson effect is not relevant under the applied intensity, see [7]. As expected for these simple cases, the agreement between analytical and FEM distributions is very good. Fig. 5 shows the vertical thermal flux vs. length for cases I (left) and coincident II, III (right). Analytical and FE results differ slightly at the edges, since the mesh is not highly refined to capture the Newmann boundary condition. For short times, the flux distribution shapes are very similar in the three cases, since the T distributions have the also the same shape. Near steady-state, a constant distribution is obtained for case I (from the quasi-linear T distribution) and a linear for case II (from a quadratic T). In addition, for case II the heat flux changes sign due to the second term in the right-side of (11)-bottom prevailing over the first one. Note that the prescribed electric flux is negative, according to the flux direction into the p-type TE. 4.2. Case IV Case IV is a non-linear and non-homogeneous hyperbolic problem. The analytical solution for the temperature is given in [8] and is based on a Laplace transform solution technique. However, several differences between this solution and our numerical results are found, thus, a new corrected comparison solution is deduced and summarized in the Appendix. Furthermore, in the present work the heat flux is calculated using a semi-analytical procedure that combines the analytical temperature solution with finite differences to obtain the gradient in (11)-bottom. No analytical results for V are generated due to the necessity of using numerical techniques in all steps. The response in this case is going to be hyperbolic at the beginning (transient state) and parabolic after some transitory time (permanent state). The first behavior is evident at early times where the V, T distributions appear as piecewise functions with sharp discontinuities near the edges. This occurs since both energies travel as finite velocity waves, see Fig. 6. In addition, strong changes of T occur when discontinuous initial boundary conditions are used, inducing a Gibbs phenomena. For larger times, the parabolic behavior dominates and, therefore, the numerical algorithm must be:  robust and efficient to automatically solve both parabolic and hyperbolic problems depending on the instant;  able to have controllable numerical dissipation to mitigate the Gibbs phenomena.

Fig. 3. p-type thermoelement geometry, applied boundary conditions, dimensions fluxes and coordinates.

For a linear problem and according to [13], no numerical oscillations will appear if the dynamic algorithm is regularized by the relationship Dt 6 Ch/v, where C is the Courant number [26] and

98

R. Palma et al. / Comput. Methods Appl. Mech. Engrg. 213–216 (2012) 93–103

Fig. 4. Voltage (left) and temperature (right) distributions along thermoelement length for cases I (top row) and coincident II, III (bottom row) at three times in [s], see Table 2. Analytical results with lines, finite element with circles.

Fig. 5. Thermal fluxes for cases I (left) and coincident II, III (right) at three times in [s], see Table 2. Analytical results with lines, finite element with circles.

pffiffiffiffiffiffiffiffiffiffiffiffiffiffi

v ¼ j=qcs,

the linear wave velocity. The cases considered in the present work are highly non-linear, therefore the spectrum of the algorithm is not enough evident for a good estimation of the ratio Dt/h. In our analyses a structured mesh of 200 elements (h = 8.1429  106 [m]) in the x3 direction is used. Since a precise spectral analysis is not an objective of this paper, a value C = 1/6 is chosen, selection based on a series of numerical tests. If smaller values of C (e.g., 1/10 and consequently smaller Dt) are used the oscillation appears again, as seen in Fig. 6 bottom. Larger values of C will artificially smear the distribution at the front itself. In order to evaluate the performance of the time integration scheme, N-b and HHT algorithms are compared using parameters  , see [27] and Table 3. that produce different spectral radii q  should stay close to unity for small to intermeAccording to [28], q b ! 1,   0:5 only when Dt= T diate time steps and decrease to q b is the undamped natural period. where T Fig. 6 top shows T vs. length x3 for an assumed s = 0.02 [s] (as in most of this section) at t = 0.06 [s]; for shorter times the comparison analytical solution does not converge since it involves Riemman inverses; the numerical one does although with increasing Gibbs noise. The analytical solution (solid line) and the FE (dashed) are obtained with N-b and HHT algorithms (these are unconditionally stable for the parameters tested when applied to linear problems).

For N-b, the use of standard parameters (test a) results in numerical oscillations since this algorithm is non-dissipative. The HHT algorithm (test c) slightly oscillates due to its low numerical damping. Therefore, in the remainder of the present work the parameters of test b will be used, since according to [27] this choice results in the highest numerical dissipation. Fig. 7 shows the V (left) and T (right) distributions along the TE at several instants. In the first, it can be observed the transition from very small times and hyperbolic solution to higher ones and parabolic, the last with almost linear distribution between boundary values. Notice in the curve for t = 0.06 [s] the electric wave fronts from left and right at around 0.28 and 0.8 [mm], with a non-linear distribution inside the wave and linear between the fronts. For t = 0.13 [s] both waves almost collide at around 0.63 [mm], point shifted to the right. This is so since the T boundary condition is higher at the left: the energy and consequently velocity are also higher for the left wave. In the T distribution a similar transition from hyperbolic to parabolic can be seen. The collision between the thermal waves is again clear for t = 0.13 [s] at the same TE point, a logical result since the thermal and electric waves are coupled. For t = 0.06 [s] it can be appreciated that inside the wave zone the slope of the distribution is higher than the obtained with a parabolic model, with the consequence of a higher energy confined in a smaller volume. The agreement between analytical and numerical results is very good,

R. Palma et al. / Comput. Methods Appl. Mech. Engrg. 213–216 (2012) 93–103

99

Fig. 8. Semi-analytical and numerical comparison of front wave thermal flux vs. thermoelement length at t = 0.06 [s] for s = 0.02 [s].

Fig. 6. Analytical temperature solution for t = 0.06, s = 0.02 [s], C = 1/6 (top). Newmark-b and different HHT parameters to study Gibbs phenomena (middle), see Table 3. Adjusting by trial and error the Courant number in Test b to avoid oscillations (bottom).

Table 3 Newmark-b and HHT parameters used for numerical testing. Test

Algorithm

 q

 b

c

a

a

b

N-b N-b N-b N-b

1 3/2 1/2 0

1/4 9/14 4/9 1

1/2 16/49 5/6 3/2

– – – –

c

HHT HHT

3/2 1/2

9/14 4/9

16/49 5/6

6/7 2/3

except for a small region around the wave front at t = 0.06 [s]; this is due to the intrinsic oscillations from the discrete time integration scheme and is very difficult to avoid. Fig. 8 compares the semi-analytical and FEM thermal fluxes for t = 0.06 [s]. Two peaks due to the discontinuous boundary conditions at the wave front are evident, representing the propagation

of this wave. Even if a nil flux at the TE center could be expected, it has a constant negative value. The reason is found in (11) bottom: the first right-term is zero due to a constant T distribution where the wave is not present (see Fig. 7). This distribution is due to the electric coupling, see the last term in (11) top. Also notice that the second term in (11) bottom will be negative due to the prescribed j. The agreement is again very good, even at the peak maximums, although it has to be considered that the two maximums are affected by the respective choices of the numerical parameters. As mentioned before, the behavior is hyperbolic in the initial instants and becomes parabolic later. This is clearly seen in Fig. 9 top, where the heat flux is plotted for several times. For t < 0.06 [s] the peaks described in the previous figure are visible, but for larger times they become smoother. At t = 0.13 [s] both waves collide, and after that the flux will be linear as in Fig. 5. The mathematical explanation of this attenuation can be seen in (11)-top: the first term in the left hand side is hyperbolic, the second parabolic and the last in the right dispersive. In the initial instants the hyperbolic term is dominant; physically this can be interpreted as a ballistic motion of either electrons or holes. The wave is constantly attenuated by the dissipation introduced by the last term up to a time in which the influence of the s terms vanish, and the first time derivative is the dominant term. The last attenuation is obviously affected by the value of the constant s, see Fig. 9 bottom, where it can be seen that for larger s the behavior will be more hyperbolic (more difficult to simulate) and vice versa. This figure is a test for the robustness of the finite element performance, but even for s = 0.06 [s] the simulation is satisfactory as long as the inequality Dt 6 Ch/v is used. 5. A complete three-dimensional simulation There are in the literature several studies aimed to the optimization of TE geometry, see [29,30] for example. However, these studies use 1D models (at most 2D), not being able to observe 3D effects. The aim of this section is to study the 3D temperature propagation inside a TE by means of the hyperbolic model validated in the previous section. Note that the propagation of V is proportional to that of T through the thermoelectric coupling, and therefore will not be included here.

Fig. 7. Voltage (left) and temperature (right) vs. thermoelement length at three times in [s]. For voltage, only finite element results; for temperature, analytical results with lines and finite element with circles.

100

R. Palma et al. / Comput. Methods Appl. Mech. Engrg. 213–216 (2012) 93–103

Fig. 9. Front wave thermal flux vs. thermoelement length at several times for s = 0.02 (top) and relaxation time influence in the wave shape (bottom) at t = 0.06, both in [s]. Only finite element results represented.

The function of a cooling thermoelectric device is to transport heat from the hot to the cold faces. As mentioned before, inside the device there are three heat fluxes (Peltier, Joule, Fourier), the first being the most interesting from a performance standpoint; Joule and Fourier fluxes are bulk effects while Peltier is a surface one. For this reason, Hoyos et al. [29] recommends to optimize the TE geometry widening the hot face section, to increase the Peltier flux and thus the device performance. To study the validity of this solution, this section analyzes a truncated pyramid TE, see Fig. 10. Sections at hot and cold faces are 1.96 and 1.58 [mm2] respectively, the length L = 0.28 [mm] and the prescribed vertical flux j = 2.65  106 [A/m2]. The temperatures at the hot and cold faces are fixed to Th = 50, Tc = 30 [°C], respectively. Initially (t = 0), the temperature is T0 = 0 °C everywhere. From an FE point of view, the element sizes and time steps are the ones optimized in the previous section. Fig. 11 shows the T propagation using the parabolic (left, s = 0) and hyperbolic (right, s = 0.02 [s]) models for four representative instants. The T propagation velocity is larger from the hot face than that from the cold, since in the former the prescribed T and the area are both larger and therefore the applied energy is higher. There are significant differences between both models. For parabolic, a T diffusion from hot to cold faces is observed that could easily be simulated by a 1D model, except for small, non-relevant

Fig. 10. Truncated pyramid thermoelement dimensions for three-dimensional simulation. On hot face surface 1.96 [mm2] and Th = 50 °C; on cold 1.58, Tc = 30 with constant prescribed vertical flux j = 2.65  106 [A/m2].

variations at the edges. For hyperbolic, many non-uniformities at several sections are observed, requiring a 3D model. The main differences between both models are:  At initial t = 7.74  103 [s]. For parabolic, the velocity of propagation is infinite, instantly increasing T by 3.47°C in all the domain except the faces. In addition, the vertical distribution is linear and constant at every section. For hyperbolic, the wavefront increases T only in the wave front, while the rest of the domain stays at the initial T0. Furthermore, T at the inclined edges increases due to wave reflection at these edges, reaching a value of 53.3 °C. This increase will reduce the performance of the truncated geometry TE.  At medium t = 3.44  102 [s]. For parabolic, the lower half increases its T and the rest reaches 25 °C all due to conduction. The vertical distribution is still constant at every section, except near the edges where concentrations occur. For hyperbolic, the wavefront from the hot face progresses upwards to more than half the TE, while the wavefront from the cold face does only one quarter. At this instant, the edge overheating is progressing with the front, and it also extends across the bottom of the TE due to the constant supply of energy. In the rest of the front section the T decreases to fulfill the conservation of total energy, see (2).  At medium t = 6.51  102 [s]. For parabolic, the system has clearly reached the steady-state with an almost linear distribution between Th and Tc. For hyperbolic, both wavefronts have already met. The concentration at the edges is now more evident, due to the multiple bounces that occur in the inclined edge.  For larger times and hyperbolic model, the wavefront reflected at the cold face returns to the hot one, reaching the steady-state T distribution, that as expected is very similar to that predicted by the parabolic. Although the geometry of this example is meso, these results could be relevant in the design and optimization of thermoelectric

101

R. Palma et al. / Comput. Methods Appl. Mech. Engrg. 213–216 (2012) 93–103

PARABOLIC MODEL

Time = 7.74E-03

Time = 3.44E-02

Time = 6.51E-02

Time = 1.29E-01

3.47E+00

HYPERBOLIC MODEL

-8.43E-09

2.30E+01

2.30E+01

2.90E+01

2.90E+01

3.50E+01

3.50E+01

4.10E+01

4.10E+01

4.70E+01

4.70E+01

5.00E+01

5.33E+01

2.78E+01

5.93E-01

2.30E+01

2.30E+01

2.90E+01

2.90E+01

3.50E+01

3.50E+01

4.10E+01

4.10E+01

4.70E+01

4.70E+01

5.00E+01

5.28E+01

3.00E+01

3.00E+01

2.30E+01

2.30E+01

2.90E+01

2.90E+01

3.50E+01

3.50E+01

4.10E+01

4.10E+01

4.70E+01

4.70E+01

5.00E+01

5.34E+01

3.00E+01

3.00E+01

2.30E+01

2.30E+01

2.90E+01

2.90E+01

3.50E+01

3.50E+01

4.10E+01

4.10E+01

4.70E+01

4.70E+01

5.00E+01

5.00E+01

Fig. 11. Temperature distribution inside the truncated pyramid thermoelement for the parabolic (left) and hyperbolic (right) models at several instants. For parabolic model s = 0, for hyperbolic s = 0.02 [s] assumed. Boundary conditions and geometry shown in Fig. 10.

Fig. 12. Analytical temperature distribution at t = 0.06 [s] for cases I–IV vs. thermoelement length. Boundary conditions and prescribed values displayed in Fig. 3.

micro-devices under fast operation modes. For instance, an electron–phonon thermal mismatch due to different thermal and electrical carriers was reported in [31] for cooling microthermoelectrics. This non-equilibrium has until now been modeled by molecular dynamics or Monte Carlo techniques; however, the present combined parabolic–hyperbolic model could simulate these effects through the introduction of an empirical value of s.

6. Conclusions This article has presented a non-linear and transient finite element formulation to simulate thermoelectric behavior under parabolic and hyperbolic heat conduction. The formulation has been implemented into the computational code FEAP. Numerically, non-linearities and transient hyperbolicity have been addressed

102

R. Palma et al. / Comput. Methods Appl. Mech. Engrg. 213–216 (2012) 93–103

by Newton–Raphson and by Newmark-b and HHT algorithms, respectively. Analytically, four 1D solutions have been developed to validate the finite element results; the first three solve parabolic problems and both numerical and analytical results agree very well, using standard parameters of Newmark-b. The fourth validation consists on a hyperbolic problem with discontinuous initial boundary conditions. In this example, V and T distributions present sharp discontinuities, and the numerical results numerical oscillations. To mitigate this Gibbs phenomena three numerical test have been tested: Newmark-b with standard parameters, HHT and New ¼ 1, c  ¼ 3=2. The latter provides the best results mark-b with b since this parametrization gives the highest numerical dissipation. In addition, the time steps and element sizes have been regularized using a linear relationship with a Courant number C = 1/6. This C has been chosen by trial and error, since the problem is highly non-linear and a spectral analysis is not one of the current work objectives. Finally a complete three-dimensional example is shown, with conclusions related with the relevance of thermoelectric pellet shapes under wave propagation. Physically, the obtained results show a T wave propagation, of thermal flux and of V (due to the thermoelectric coupling). This numerical tool could be applied to the design of micro-devices for novel applications with fast processes, where the second sound effect could be relevant. Acknowledgments This research was partially supported by the Spanish Ministry of Education through Grant No. FPU AP-2006-02372 and also from Grants MICINN BIA-2008-00522, CSD2008-00037 Canfranc Underground Physics, Excelencia Junta Andalucía P08-TEP-03641 and ‘‘Ayudas Investigación’’ from UPV. The authors would also like to thank Prof. Guillermo Rus for his valuable contributions. Appendix A. Analytical solutions The analytical solution of case III is directly given in [25]. In this appendix specially developed solutions for I, II, IV are presented. The T distribution along the TE for I, II are reported in [24]. The V solution is calculated with (9), the thermal flux one by (11)bottom. Case I

  1 P Vðx3 ; tÞ ¼ a DLT x3  DT þ cn SxeKt ; n¼1

1 P

DT L

Tðx3 ; tÞ ¼ T c þ x3 þ cn SxeKt ; n¼1   1 P q3 ðx3 ; tÞ ¼ j DLT þ cn Cx nLp eKt ;

ðA:1Þ

n¼1

where DT = Th  Tc, K = b(np/L)2, b = j/qc, Sx  sin(npx3/L), Cx  cos(npx3/L) and L is the TE length. In addition

2 cn ¼ L

Z

L

0



 DT x3 Sxdx3 : Tc  L

ðA:2Þ

Case II

Vðx3 ; tÞ ¼ jc3 Lx þ a A2 L

 DT L

 1

P  A1 x3 Lx cn SxeKt ; n¼1 1 P

x3 þ T C  þ cn SxeKt ; n¼1   1 P q3 ðx3 ; tÞ ¼ j AL2  2A3 x3 þ cn Cx nLp eKt þ aTðx3 ; tÞj3 ;

Tðx3 ; tÞ ¼

A3 x23

n¼1

ðA:3Þ

2

2

where Lx = L  x3, A1 = j3/2jc, A2 ¼ DT þ j3 L2 =2bqcc, A3 ¼ j3 =2bqcc and

cn ¼

2 L

Z

L

0

  A2 x3 þ T c  A3 x23 Sx dx3 : L

ðA:4Þ

Case IV The T distribution along the TE was reported in [8] using the Laplace transform technique. However, we believe there are several errors in this work, for which the result is developed again. From the reference, we use the dimensionless parameters n = x3/L and g = tj/ qcL2 but (10) is rewritten in an amenable form to be solved by the Laplace technique

@ nn fTg þ F 1 @ n fTg þ F 2  @ g fTg  sF 3 @ gg fTg n o þsF 4 @ g F 2 ALj2  tajLA2 F 1 @ n fTg ¼ 0;

ðA:5Þ

h

Tðn; 0Þ ¼ 0;

Tð0; gÞ ¼ T h ;

Tð1; gÞ ¼ T c :

Note that the boundary conditions have been included. Regarding [8], F1 = thIL/jA, F2 = I2L2/ jcA2 have been modified and F3 = j/qc L2, F4 = 1/ qcA are added. Applying the Laplace transform

Wðn; sÞ ¼ P11 ek11 n þ P12 ek12 n þ sðsF Fs22 þsÞ ; 3

WðnÞ ¼ 0;

Wð1Þ ¼ Tsh ;

Wð0Þ ¼ Tsc ;

ðA:6Þ

where s is time in the Laplace domain. The boundary conditions in (A.6) are different from those in [8]: in this reference they are not transformed. Then, the following parameters are modified Th k12 s þe

P11 ¼



F2 T c sðsF 3 s2 þsÞ s





s

F2

ðsF3 s2 þsÞ

ek11 ek12

;

P12 ¼ Tsc  P11  sðsF Fs22 þsÞ ; 3

 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi k11 ¼ 12 A1 A21  4B1 ; k12 A1 ¼ F 1 ðt h L2 þ AajssF 4 Þ;

ðA:7Þ

B1 ¼ t 2h L4 ðs þ F 3 s2 sÞ:

Finally, (A.6) are inverted using the Riemann-sum approximation as in [8]. In order to compare the different analytical solutions, Fig. 12 shows the T distributions obtained at t = 0.06 [s], using the geometry and boundary conditions shown in Fig. 3. Only T distributions are shown since no analytical solutions for V and heat flux distributions for case IV are available. Results for cases II, III are very similar, since the Thomson effect is not relevant under the applied electric intensity. For cases I–II, the distributions are quadratic due to the boundary conditions. The small T increase for cases II, III with respect to I is due to the Joule effect, which acts as a heat source. For case IV, the T propagates as a wave, confining all energy into the wavefront; the T wave has not reached the rest of the TE and remains at the initial 0 °C. References [1] S. Riffat, X. Ma, Thermoelectrics: a review of present and potential applications, Appl. Therm. Engrg. 23 (2003) 913–935. [2] G. Snyder, J. Fleurial, T. Caillat, Supercooling of Peltier cooler using a current pulse, J. Appl. Phys. 92 (3) (2002) 1564–1569. [3] C. Cattaneo, Sulla conduzione del calore, Atti Seminario Mat. Fis. Univ. Modena 3 (1938) 83–101. [4] D. Jou, G. Lebon, Extended Irreversible Thermodynamics, Springer-Verlag, Berlin, Heidelberg, 1996. [5] G. Lebon, D. Jou, J. Casas-Vázquez, Understanding Non-equilibrium Thermodynamics, Springer-Verlag, 2008. [6] E. Antonova, D. Looman, Finite elements for thermoelectric device analysis in ANSYS, in: International Conference on Thermoelectrics, 2005. [7] J. Pérez-Aparicio, R. Taylor, D. Gavela, Finite element analysis of nonlinear fully coupled thermoelectric materials, Comput. Mech. 40 (2007) 35–45.

R. Palma et al. / Comput. Methods Appl. Mech. Engrg. 213–216 (2012) 93–103 [8] M. Alata, M. Al-Nimr, M. Naji, Transient behaviour of a thermoelectric device under the hyperbolic heat conduction model, Int. J. Thermophys. 24 (6) (2003) 1753–1768. [9] J. Gembarovic, V. Majernik, Non-fourier propagation of heat pulses in finite medium, Int. J. Heat Mass Transfer 31 (5) (1988) 1073–1080. [10] H. Chen, J. Lin, Numerical analysis for hyperbolic heat conduction, Int. J. Heat Mass Transfer 36 (11) (1993) 2891–2898. [11] A. Barletta, E. Zanchini, Hyperbolic heat conduction and local equilibrium: a second law analysis, Int. J. Heat Mass Transfer 40 (5) (1997) 1007–1016. [12] M. Manzari, M. Manzari, On numerical solution of hyperbolic heat conduction, Commun. Numer. Methods Engrg. 15 (1999) 853–866. [13] R. Ciegis, Numerical solution of hyperbolic heat conduction equation, Math. Modell. Anal. 14 (1) (2009) 11–24. [14] H. Gómez, I. Colominas, F. Navarrina, M. Casteleiro, A discontinuous Galerkin method for a hyperbolic model for convection–diffusion problems in CFD, Int. J. Numer. Methods Engrg. 71 (2007) 1342–1364. [15] S. Bargmann, P. Steinmann, Theoretical and computational aspects of nonclassical thermoelasticity, Comput. Methods Appl. Mech. Engrg. 196 (2006) 516–527. [16] R. Taylor, FEAP, A Finite Element Analysis Program: User Manual, University of California, Berkeley, 2010. . [17] H. Hilber, T. Hughes, R. Taylor, Improved numerical dissipation for time integration algorithms in structural dynamics, Earthquake Engrg. Struct. Dyn. 5 (1977) 282–292. [18] J. Jackson, Classical Electrodynamics, John Wiley and Sons Inc., 1962. [19] A. Eringen, G. Maugin, Electrodynamics of Continua I, Springer-Verlag New York, Inc., 1990.

103

[20] A. Chakraborty, K. Ng, Thermodynamic formulation of temperature-entropy diagram for the transient operation of a pulsed thermoelectric cooler, Int. J. Heat Mass Transfer 49 (2006) 1845–1850. [21] J.L. Pérez-Aparicio, R. Palma, R.L. Taylor, Finite element analysis and material sensitivity of Peltier thermoelectric cells coolers, Int. J. Heat Mass Transfer, 2011, doi:10.1016/j.ijheatmasstransfer.2011.08.031. [22] O. Zienkiewicz, R. Taylor, J. Zhu, The Finite Element Method: The Basis, Elsevier/Butterworth-Heinemann, 2005. [23] MELCOR, Application notes for thermoelectric devices, Tech. Rep., Laird Technologies, Melcor. . [24] A. Polyanin, Handbook of Linear Partial Differential Equations for Engineers and Scientist, Chapman & Hall CRC, 2002. [25] W. Seiffert, M. Ueltzen, E. Muller, One-dimensional modelling of thermoelectric cooling, Phys. Stat. Sol. A 194 (1) (2002) 277–290. [26] K.F.R. Courant, H. Lewy, Uber die partiellen differenzengleichungen der mathematischen physik, Math. Annal. 100 (2007) 32–74. [27] T. Hughes, The Finite Element Method, Linear Static and Dynamic Finite Element Analysis, Prentice-Hall, Inc., 1987. [28] T. Fung, Numerical dissipation in time-step integration algorithms for structural dynamic analysis, Prog. Struct. Engrg. Mater. 5 (2003) 167–180. [29] G. Hoyos, K. Rao, D. Jerger, Fast transient response of novel Peltier junctions, Energy Convers. 17 (1977) 23–29. [30] Y. Cheng, W. Lin, Geometric optimization of thermoelectric coolers in a confined volume using genetic algorithms, Appl. Therm. Engrg. 25 (2005) 2983–2997. [31] L. da Silva, M. Kaviany, Micro-thermoelectric cooler: interfacial effects on thermal and electrical transport, Int. J. Heat Mass Transfer 47 (2004) 2417– 2435.