Device Efficiency of TEG and TEC :
      One-dimensional  Modelling of
      Thermoelectric Cooler and Generator

Dr. W. Seifert, S. Walczak, Dept. of Physics, University Halle-Wittenberg
           and
          Dr. E. Mueller,  Institut  f"ur Werkstoff-Forschung, DLR, Cologne, Germany

Introduction and theoretical background

Thermoelectric effects are discussed here within the framework of the linear Onsager 
               theory. Under isotropic conditions, the constitutive relations are
                                   (Overscript[jel, →]) = σ (Overscript[E, →]) - σ S  ∇T    ;     (Overscript[jq, →]) = S T  (Overscript[jel, →])  - κ  ∇T                                 (1)
               with the electric field (Overscript[E, →]), electric current density (Overscript[jel, →]) ,  temperature T,  heat flux  (Overscript[jq, →])  ,
               electric conductivity σ,  thermal conductivity κ ,  Seebeck coefficient S and Nabla ∇.
               Assuming steady-state conditions, the gouverning equations result  from the
               principles of conservation of charge and energy
                    
                   ∂U/∂t+∇.(Overscript[((Subscript[j, q])), ⇀])=(Overscript[((Subscript[j, el])), ⇀]) .  (Overscript[E, ⇀])      =>   steady-state case:        ∇.(Overscript[((Subscript[j, q])), ⇀])= (Overscript[((Subscript[j, el])), ⇀]) .  (Overscript[E, ⇀])       (2)
                   ∂ρ / ∂t + ∇.(Overscript[((Subscript[j, el])), ⇀])  =  
        0             =>   steady-state case :      ∇.(Overscript[((Subscript[j, el])), ⇀])  =  
        0  ,            
                    
               The coupling between heat transport and charge transport is caused by two 
               (inverse) thermoelectric effects: Heat is transported by a flowing current or is used
               to generate electricity based on the Peltier and Seebeck effects, respectively. 
               
               An analytical solution of the transport of charge and heat through a thermoelectric 
               element can be found for  the 1D case of constant material properties (1D-CPM).
               Concerning Peltier cooling, we refer to our papers :
                     W. Seifert, M. Ueltzen, E. Mueller: physica status solidi(a) 194, 1, 277-290 (2002)
                     E. Mueller, S.Walczak, W. Seifert: physica status solidi(a) 203, 8, 2128-2141 (2006)
                     S.Walczak, W. Seifert, E. Mueller: Mater. Sci. Forum Vols. 492-493, 
                           Trans Tech Publications,  Switzerland (2005) pp. 507-515.
                     E. Mueller et. al.: Proc. 24th Int. Conf. on Thermoelectrics, Clemson, SC, USA, 
                          June 19-23, 2005, IEEE, Piscataway, NJ (2006) pp. 352-357.
                         
                A central target of theoretical TE studies is to elaborate recipes for optimum design 
                of TE elements integrated into TE modules. However, clear guidelines have not been 
                available yet for long.
                The new compatibility concept introduced for thermoelectric generators by
                     G. J. Snyder, T. S. Ursell,  Phys. Rev. Letters 91, 14, 148301- 1-4 (2003)
                provides a new approach based on first principles of thermoelectricity and
                offers a new pathway for a rational search of efficiently graded TE systems.
               
               In this notebook, this concept is applied for constant material parameters  S, κ, σ.
               Closely following the approach suggested by Snyder and Ursell, thermoelectric generator 
               (TEG) and cooler (TEC) will be treated here in the same mathematical formalism. 
               In both cases, a 1D model is used considering a single p-type pellet with a constant 
               cross sectional area. Furthermore, the electric current density jel  is assumed to be a 
               constant. 
               How to generalize the compatibility concept also to Peltier cooling, is shown in
                    W. Seifert, E. M"uller, S.Walczak:  Generalized analytic description of one-dimensional
                    non-homogeneous TE cooler and generator elements based on the compatibility approach,
                    Proc.  25th Int. Conf. on Thermoelectrics (ICT 2006), Wien, AU, 6.-10.8.2006. 
               
               For clarity we will refer to the heat absorbing side (cooling side, mostly cold side of the 
               cooler and hot side of the TEG, respectively) and the sink side. TEG and TEC cases are
               then mainly distinguished by opposite sign of the temperature difference !
                
               Without any doubt, global maximisation of a performance parameter (electric output power 
               or efficiency of a TEG; cooling power or coefficient of performance (C.O.P.) of a Peltier device) 
               is a suitable guideline for an empirical approach to TE device optimisation. In this notebook,
               these parameters are computed within the frame of the Constant Properties Model (1D-CPM).
               Besides, the relative current density  u(T) and the compatibility factor  s(T) are plotted.

Definitions used below : 

               Generator:    Efficiency  η = net electrical power output / thermal power supplied = -Pelek / H  
              
                   Cooler:     Coefficient of performance C.O.P. = cold side cooling power  / electrical power supplied = H / Pelek

Input data

Clear["Global`*"] 
 Off[General::"spell1"];

Input  parameters  ( SI-units ! ) :  

                                                       L :   length of the active material
                                             T1 ≡ Ta :  absorbing side temperature  T(x=0) = Ta  in Kelvin 
                                             T2 ≡ Ts :  sink side temperature           T(x=L) = Ts  in Kelvin             
                                             
                                             T2 ≠ T1   ⟹  Generator :  T1 > T2  ;   Cooler :  T2 > T1   !

L =  0.005         
        T1  = 400             (* 
      T1  =  Ta ⟹ 
            absorbing side *)     
                                      (* 
      alternative test,  
      for Cooler :  T1  =  260  *) 
        T2  = 300             

  matparam = { 
        S -> 180*(10^(-6)),  κ -> 1.35, σ -> 
          140000};     (((S^2) σ)/κ) /.  
      matparam ;
  matparam=Append[matparam, 
        Z -> %]

Calculate Z·T  at  sink side temperature T2

Z*T2  /.  matparam

General hint:
                    Some of the graphics below are optimized to the set of parameters used.
                    For variations of parameters use another  PlotRange or  "PlotRange->All"  if necessary !

Loop over current density jel  : predefinitions

Nmax = 
        20;                 

    If[T1  >  T2,  
        jopt = ((S*σ*(T1 - T2))/(2*L))   /.  matparam , 
        dummy]  ;  
    If[T1  <  T2,  
        jopt = ((T2*S*σ)/L)   /.  matparam ,  
        dummy]  ;
    If[T1  >  T2,  
        Print[" TEG: recommented value for jelstep ≤ ",  
          jopt/10], dummy] ;
    If[T1  ≤   T2, 
        Print[" TEC: recommented value for jelstep ≤ ",  
          jopt/20],  
        dummy] ;      

   jelstep = 
        25000  ;                     

    Print[" input value for jelstep = ",  
      jelstep]

Pay attention to the recommented value of jelstep giving you an orientation for TEG / TEC !

Constant Properties Model for  TEG / TEC

Print[" Loop plots "];
 Print["             first plot:  black: temperature profil  T(x) \
"];
\(\ \(Print["             second plot:  blue: u(T) ; red: s(T) "];\)
\)
 xy = 
        Array[dat, {Nmax, 
            4}];

\(\ \(For[\ inx = 1, \ 
        inx < \((Nmax + 
              1)\), \ \(inx++\), \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \  (*\ \ \ Loop\ over\ current\ density\ jel\ \ *) 
\ \ \ \ \ \ \ jel = jelstep*
            inx\ \ ; 
\ \ \ \ \ \ \ \ \ \ xy[\([inx, 
              1]\)]\  = \ 
          jel\ ; 
co = \(\(jel\^2\)\(\ \)\)\/\(σ*κ\)\  /. \ matparam; 
\ \ \ \ dT\  = \ 
          T2 - T1; 

\ \ \ \ \ \ \ \ \ \ \ c1 = dT/L + \ co/2*L\  /. \ 
            matparam; \ \ \ \ \  (*\ \ c1\  = \ \(gradT\ at\ x = 
              0\)\ \ *) 

\ \ T[
            x_] = \(-co\)/2*x^2 + c1*x + 
            T1; 

\ \ \ \ If[T1 > T2, \ 
            Tstart = T2; Tend = T1, \ Tstart = T1 - 10; 
            Tend = T2 + 
                100, \ \ ]

\ \ \ \ \ Plot[T[x], {x, 0, L}, 
            PlotRange -> All\  (*\ 
              PlotRange -> {Tstart, 
                  Tend}\ *) \ ]; 

\ \ \ \ \ \ \ \ \ \ uc = \(jel/κ\)/
            c1; \ \ \ \ \ \ \ \ \ 

\ \ \ \ \ \ \ \ \ u[T_] = 
          1/\@\(1/uc^2 - 2*κ/σ*\((T - T1)\)\)\ \  /. \ 
            matparam\ ; \ \ \ \ 

If[
          T1 > T2, \ \ \ \ \ 
\ \ \ \ \ \ \ \ \ \ sTEG[
              T_] = \(\@\(1 + \ T*Z\) - 1\)\/\(S*T\)\  /. \ 
              matparam; 
\ \ \ \ \ \ \ \ \ Plot[{u[T], 
              sTEG[T]}, {T, T2, T1}, 
            PlotRange -> {0, 
                30}, 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ PlotStyle -> {RGBColor[0, 0, 1], \ 
                RGBColor[1, 0, 
                  0]}]\ \ , \ 

\ \ \ \ \ \ \ \ \ sTEC[T_] = \(\(-1\) - \@\(1 + \ T*Z\)\)\/\(S*T\)\  /. \ 
              matparam; 
\ \ \ \ \ \ \ Plot[{\(-u[T]\), 
              sTEC[T]}, {T, T2, T1}, 
            PlotRange -> {\(-100\), 
                0}, 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ PlotStyle -> {RGBColor[0, 0, 1], \ 
                RGBColor[1, 0, 
                  0]}]\ , \ \ ]; 

jq[
            x_] = \(-κ\)*D[T[x], x] + S*jel*T[x] /. 
            matparam; 
\ 
        H = jq[x]\  /. \ x -> 0; 
\ 
        PelekTE = 
          jel*S*dT /. \ matparam\ \ \  // \ 
            N\ \ ; \ \ \ \ \ \ \ \ \ \ \ \ \ \ \  (*\ \ TEG : \ 
            dT\  < \ \(0\ !\)\ \ *) 
\ 
        PelekV = 
          jel^2/σ*L /. \ 
            matparam\ ; 
\ \ \ \ \ \ \ \ \ Pelek = 
          PelekTE + PelekV; 

If[\ 
          T1 > T2, \ \ 
\ \ \  (*\ \ \ Generator\ \ \ *) 
\ \ \ \ \ \ \ \ \ \ \ η = \ \(-Pelek\)/
              H; 
\ \ \ \ \ \ \ \ \ \ \ ηglobal = \ 
            jel*\ U/H\ \ \  /. \ 
              matparam, \ \ \ \ \ \ \ \  (*\ \ global\ Def . \ 
                of\ η\ \ *) 
\ \ \  (*\ \ \ Cooler\ \ \ \ \ *) \ \ \ 
\ \ \ \ \ \ \ \ \ \ \ η = \ 
            H/Pelek; 
\ \ \ \ \ \ \ \ \ \ \ ηglobal = \ H/\((jel*
                  U\ )\)\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \  (*\ \ global\ Def . \ of\ C . O . 
                P . \ \ *) 
\ \ ]; 

\ \ \ \ \ \ \ xy[\([inx, 2]\)] = \ 
          H; 
\ \ \ \ \ \ \ xy[\([inx, 
              3]\)] = \ \(-Pelek\)\ ; \ \ \ \ \ \ \ \ \ \ 
\ \ \ \ \ \ \ \(xy[\([inx, 
                4]\)] = \ \ η\ ;\)\ \ \ \ \ \ \ \ \ \ \ 

 (*\ 
          Print["Outputs for control:"]; 
\ \ Print["  jel =", \ jel]; 
\ \ Print[" H =", \ 
            H, "  U*jel =", \ U*jel, \ " Pelek =", 
            Pelek, "  η =", \ η\ ]\ ;\ \ \ \ \ *) 
]\ ;\)
\)
 Print[" Results " ];

 Print[" jel          H            -Pelek         η / C.O.P.  \
   "] ;
TableForm[xy]

Results

Interpolation : search for maximum

tabtmp = Array[dat, {Nmax, 2}]; 
\(\(Do\ [\ \ \ \(tabtmp[\([i, 1]\)]\  = \ 
            xy[\([i, 1]\)]/
              10^6;\)\ \ \ \ ; \ \ \ 
	\ \ \ \ \ \ \ \ \ tabtmp[\([i, 
              2]\)]\  = xy[\([i, 4]\)], \ {i, 1, 
          Nmax}\ ]\ ;\)\ 
\)
\(\(If[T1\  > \ T2, \ 
        Print[\ "  η(jel)  "], \ \ Print[\ "  COP(jel) "], \ ]\ ;\)\)
plot1 = 
        ListPlot[tabtmp, 
          PlotJoined -> True];
If[T1  >  T2,  
        Do [   tabtmp[[i, 2]]  = xy[[i, 3]],  {i, 1, 
            Nmax} ] , 
                        Do [   tabtmp[[i, 2]]  = xy[[i, 2]],  {i, 1, 
            Nmax} ] ,  ];

\(\(If[T1\  > \ T2, \ 
        Print[\ " -Pelek(jel)  "], \ \ Print[\ "  H(jel) "], \ ]\ ;\)\)
plot2 = 
        ListPlot[tabtmp, 
          PlotJoined -> True];

tabeta = 
        Array[dat, {Nmax, 2}]; 
\(\(Do\ [\ \ \ tabeta[\([i, 1]\)]\  = \ 
          xy[\([i, 1]\)]; \ \ \ 
\ \ \ \ \ \ \ \ \ \ \ tabeta[\([i, 
              2]\)]\  = xy[\([i, 4]\)], \ {i, 1, 
          Nmax}\ ]\ ;\)\ 
\)
dummywert = 0;
dummyindex = 0;
Do[
        If[dummywert ≤ tabeta[[i, 2]], 
          dummywert = tabeta[[i, 2]]; dummyindex = i], {i, 1, 
          Length[tabeta]}];

etaInterp = 
        Interpolation[tabeta, 
          InterpolationOrder -> 
            2];

ResList = 
        FindMinimum[-etaInterp[x], {x, tabeta[[dummyindex, 1]], 
            tabeta[[1, 1]], 
            tabeta[[Nmax, 
                1]]}];

etamax = -ResList[[1]] 
 jipoleta = x /. ResList [[2]]

tabeta = 
        Array[dat, {Nmax, 2}]; 
\(\(If[
        T1 > T2, 
\ \ Do\ [\ \ \ \(tabeta[\([i, 1]\)]\  = \ 
              xy[\([i, 
                  1]\)];\)\ \ \ \ ; \ \ \ 
	\ \ \ \ \ \ \ \ \ tabeta[\([i, 
                2]\)]\  = xy[\([i, 3]\)], \ {i, 1, 
            Nmax}\ ]\ , 
\ 
        Do\ [\ \ \ \(tabeta[\([i, 1]\)]\  = \ 
              xy[\([i, 
                  1]\)];\)\ \ \ \ ; \ \ \ 
	\ \ \ \ \ \ \ \ \ tabeta[\([i, 
                2]\)]\  = xy[\([i, 2]\)], \ {i, 1, 
            Nmax}\ ]\ , \ \ ];\)\ \)
     dummywert = 0;
dummyindex = 0;
Do[
        If[dummywert ≤ tabeta[[i, 2]], 
          dummywert = 
            tabeta[[i, 
                2]]; 
                dummyindex = i, ], {i, 1, 
          Length[tabeta]}];

etaInterp = 
        Interpolation[tabeta, 
          InterpolationOrder -> 
            2]; 

ResList = 
        FindMinimum[-etaInterp[x], {x, tabeta[[dummyindex, 1]], 
            tabeta[[1, 1]], 
            tabeta[[Nmax, 
                1]]}];

  Pmax = -ResList [[1]] 
 jipolP = x /. ResList[[2]]

Print["  Summary :  Results of \
Interpolation "];

\(\(If[
        T1\  > \ 
          T2, 
\ \ Print[" Generator:  Max[eta] = ", \ 
          etamax, " at jopt = ", 
          jipoleta], 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Print[" Cooler: Max[COP] = ", \ etamax, "  at jopt = ", 
          jipoleta]\ \ , \ \ ]\ \ \ ;\)\)
\(\(If[
        T1\  > \ 
          T2, 
\ \ Print[" Generator:  Max[Pelek] = ", \ Pmax, " at jopt = ", 
          jipolP], 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Print[" Cooler:   Max[H] = ", \ Pmax, "  at jopt = ", 
          jipolP]\ \ , \ ]\ \ \ ;\)\)