One-dimensional  Modelling of
     Thermoelectric cooling (I)

Dr. W. Seifert, S. Walczak,  FB Physik, Universit"at Halle-Wittenberg
         Dr. E. Mueller,  Insitut f"ur Werkstoff-Forschung,  DZfLR,  K"oln

Introduction and theoretical background

Thermoelectric effects are caused by coupling between charge transport and heat 
               transport. Under  isotropic conditions, within the frame of the linear Onsager theory, 
               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  ,            
                    
               Assuming steady-state conditions and applying an one-dimensional  formulation 
              with a constant electric current density jel through a constant cross-sectional  
               area Aq ( jel > 0 along the positive x-axes, active thermoelectric material in the
               region  0 ≤ x ≤ Lprobe ),  the basic equations of our 1D model result from Eqn. (1)
               and (2) :
                                 ∂/∂x(-κ ∂T/∂x + S T jel) = jel^2/σ + S jel  ∂T/∂x                             (3)
                                   
                                   E ≡ -  ∂(Subscript[φ, el])(x)/∂x  =  jel/σ + S  ∂T/∂x                                            (4)
               
                Assuming the simplifying case of constant material parameters  S, κ, σ, 
                Eq. (3) reduces to the simple differential equation  
                                   (T^″)(x) = - c0 ;   c0 =  jel^2/(σ*κ)                                                   (5)
               for the resulting temperature profile over the length Lprobe of the thermoelectric
               material. The solutions of Eq. (5) and Eq. (4)     
               
                                T(x) = -c0/2*x^2 + c1*x +c2 ,                                                       (6)
                              (Subscript[φ, el])(x) =  S*c0/2*x^2 - x*( jel/σ + S*c1 ) + c3                               (7)
                                              
              with  constants       c1 = 1/L* (T2 - T1 + c0/2*L^2 )   ;     c2 = T1  ;     c3 = 0,
              fulfills boundary conditions of  1. kind
                           T(x=0) = Tc  ;     T(x=L)  = Th   ;    (Subscript[φ, el])(x=0) = 0.
              Starting from Eqn. (6) and (7), the notebook gives a short overview 
              on Peltier effect and thermoelectric cooling.
              For more detailed information we refer e.g. to our paper published in
                           physica status solidi(a) 194, 1, 277-290 (2002)

Clear["Global`*"] 

     Off[General::"spell1"];

Input  parameters  (SI-units) :  
                                        Lprobe :   length of the active material
                                                Th :  hot side temperature  T(x=L) = Th   (internal parameter T2)               
                                               Tc :  cold side temperature  T(x=0) = Tc  (internal parameter T1)

Lprobe = 0.001 ;

               Th = 300 ;

               Tc = 270;

Calculate  thermoelectric figure of merit   Z = S^2*σ/κ   and add to the list of parameters

param = {L -> Lprobe,  T2 -> Th, 
          S -> 210*(10^(-6)),  κ -> 1.35, σ -> 
            65000};
               ((
(S^2) σ)/κ) /. 
        param;

param = Append[param, Z -> %]

Calculate Z·T  at hot side temperature Th

Z*Th  /.  param

General hint:
           Some of the graphics below are optimized to the set of parameters used.
           For variations of parameters change the "PlotRange"  if necessary !

Constant Properties Model (CPM)

Constant Properties Model will be applied in this notebook.
           Define and solve differential equation together with boundary conditions of  1. kind

deq1 = (T^″)[
            x] == -c0;

     eqsol =   DSolve[deq1, T[x], x]  // 
          Flatten ;

     T[x_] =  T[x]  /.  eqsol

Koeffs = 
      Flatten[Solve[{T[0] == T1,  T[L] == T2}, {C[1], 
            C[2]}]]

T[x_] =  T[x]  /. Koeffs  /. 
        c0 -> ((jel^2)/(κ σ))

calculate heat flux

jq[x_] = -κ*D[T[x], x] + jel*S*T[x]  // Simplify

calculate electric potential  φel(x)  and  prove   φel(0) = 0

φel[x_] = 
      Integrate[-jel/σ  -  S*D[T[x], x], x] // 
        Simplify
φel[0]  // Simplify

Example: for a defined value of  jel=1500000 :  plot  T(x) ,  φel(x)  and   jq(x)

            note  SI-units :          jel = 1500000 A/m^2      = 1.5 A/mm^2 ,
                                                 jq =   100000 J/(m^2  s) = 100000 W/m^2 =  0.1 W/mm^2

j0 = 
        1500000 ;  

  

\(\(\(Plot[\ \ \(\(T[x]\  /. \ param\)\  /. \ T1 -> Tc\)\  /. 
          jel -> j0, 
\ \ \ \ \ \ \ \ \ \ \ \ \ {x, 0, 
          Lprobe}, 
        PlotRange -> {250, 
            310}\ , \ \ \ \ \ \ \ 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ AxesLabel -> {" x ", " T(x)/K"}]\ ;\)\(
\)
    \)\)
jqth[
          x_] = -κ*D[T[x], x] /.  param   /. 
            T1 -> Tc  /. jel -> j0;
  jqP[x_] = jel*S*T[x] /.  param   /.  T1 -> Tc  /. 
          jel -> j0;

Plot[φel[x]  /.  param /.  T1 -> Tc  /. 
          jel -> j0, {x, 0, 
          Lprobe}, 
          AxesLabel -> {" x \
", " φel(x)/ V"}] ; 

\(\(Plot[\ {\(\(jq[x]\  /. \ param\)\ \  /. \ T1 -> Tc\)\  /. 
            jel -> j0, 
\ \ \ \ \ \ \ \ \ jqth[x], 
          jqP[x]\ \ }, \ {x, 0, 
          Lprobe}, 
\ \ \ \ \ \ \ \ \ \ \ PlotStyle -> {RGBColor[0, 0, 0], 
            RGBColor[1, 0, 0], 
\ \ \ \ \ \ \ \ RGBColor[0, 
              0, 1]\ }, 
        PlotRange -> {\(-110000\), 
            160000}, \ \ \ \ \ \ \ \ \ \ \ 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ AxesLabel -> {" x ", " jq(x)/ W/m^2"}];\)\)
jQabsorbed = jq[0] /.  param /.  T1 -> Tc  /. 
          jel -> j0;
Print["absorbed heat at the cold side x=0 : ", 
        jQabsorbed];

Last graphics shows total heat flux and its components :
                 total heat flux (black line) ;  Peltier heat flux (blue line) ;   thermal heat flux (red line)

Define cold side heat flux  jq(x=0)  in dependence on the electric current jel
          and the cold side temperature T1:

Clear[jel, T1];
jq0[jel_, T1_] = jq[x] /. param   /. x -> 0

3D plot of  jq(x=0)⩾ 0 ;  the curved line  jq(x=0)= 0  represents maximum cooling !

\(\(Plot3D[
        jq0[jel, T1], {jel, 1, 3*10\^6}, {T1, 300, 
          230}, 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ PlotRange -> {0, 120000}, 
        PlotPoints -> 
          50, 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ AxesLabel -> {\*""(Subscript[j, el]) "", \*""(Subscript[T, 1])"", \*""(Subscript[j, q])    ""}];\)\)

Example:  2D-plot of  jq(x=0,Tc=270 K)

\(\(Plot[jq0[jel, 270], {jel, 15000, 6*10\^6}, 
        AxesLabel -> {\*""  (Subscript[j, el])"", \*""(Subscript[j, q])""}, 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ PlotRange -> {\(-20000\), 80000}, 
        AxesOrigin -> {0, 0}];\)\)

Arbitrary combinations of parameters are not usefull to understand thermoelectric cooling ;
           that is why we proceed with a systematic discussion of operation conditions.

Systematic investigation of maximum Peltier cooling

Strategy :
           The thermal flux absorbed at the cold side  x = 0  
                          jq(x=0) = - c1*x + jel*S*T1  ;   c1 = 1/L*(T2-T1 + c0/2*L^2)  ;   c0 =jel^2/kappa/sigma
            shall be rewritten as a function of the temperature difference  dT=T2-T1 as follows :
                          jq(x=0) = - (kappa/L + jel*S)*dT + jel*S*T2- jel^2*L/2/sigma
            A discussion of this equation will make the situation more clear ...

Calculate heat flux  jq(x=0) at the cold side

jqx0 = jq[0]  /.  T1 -> T2 - dT  // Simplify

For maximum cooling, i.e. (Subscript[j, q])(0)-> 0, an equation for  dT  holds

eq2 = Solve[jqx0 == 0, dT] // Flatten

tempdiff[jel_] = dT /. eq2  [[1]]

Temperature difference depends on the electric current.
           Search for maximum

Solve[D[tempdiff[jel], jel] == 0, jel]   // Flatten

Rewrite this solution by means of the thermoelectric figure of merit   Z = (S^2)σ / κ

eq3 = 
      FullSimplify[
          Solve[D[tempdiff[jel], jel] == 0, jel], {κ > 0, 
            Z == (((S^2) σ)/κ)}] // Flatten

As jel  flows along the positive x-axes ( jel > 0 ), only  the second term will be used

jelopti = jel /. eq3 [[2]]

Calculate maximum temperature difference   Δ(Subscript[T, max]) and  minimum  cold side temperature  (Subscript[T, 1  
                  min]):

dTmax = 
      FullSimplify[
          tempdiff[jelopti], {κ > 0, 
            Z == (((S^2) σ)/κ)}] // 
        Expand
T1min = Simplify[T2 - dTmax] // Expand

values ...

jelopti /. param  // N
 dTmax   /. param  // N
 T1min /. param  // N

Temperature profile in the case of  maximum cooling

tempmaxc[
        x_] = T[x]  /.   jel ->  jelopti   /. 
          T1 -> T1min   /.  param  

      
\(\(\(pl1 = 
        Plot[tempmaxc[x], {x, 0, Lprobe}, 
          AxesLabel -> \ {"x", "T(x)"}, 
\ \ \ \ \ \ PlotRange -> {220, 340}];\)\(
\)
    \)\)
tempmaxc[0]

Hint :
             Alternatively, in the case of maximum cooling, temperature profile and electric potential can be defined
             using the structure of the differential equation's solution ( see Eqn. (6),(7) ):

Note:
             A comparison of the expressions for  jelopti (see Out[33])  and  T1min (see Out[35]) leads to the 
             well-known relation   
                                                                               jelopti = ((σ  S)/L) T1min 
             Proof :

jelopti  - σ*S/L*T1min    /. 
        Z -> (((S^2) σ)/κ)   //Simplify

Further note, that the equation    jopt = ((sigma  S)/L) T1   is also valid  in the case of isothermal outside
             conditions  T1=T2, i.e. the optimum current is, in general, proportional to the cold side temperature !
             
            We prove this in analogy to the procedure used above:

f[jel_] = jqx0  /.  dT -> 0  // Simplify

Solve[D[f[jel], jel]  == 0, jel]   // Flatten  // 
      Simplify

Owing to the isothermal condition T2 = T1 , the result is now :   jelopti2 = ((σ S)/L) T2 .

jelopti2 = jel  /.  %

Calculate and plot temperature profile in the case of isothermal pumping

jelopti2  /.  param
tempiso[
        x_] = T[x]  /.   jel ->  jelopti2   /. 
          T1 -> T2   /.  param  

      
\(\(pl2 = 
        Plot[tempiso[x], {x, 0, Lprobe}, 
          AxesLabel -> \ {"x", "T(x)"}, 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ PlotStyle -> Dashing[{0.02}], 
          PlotRange -> {220, 340}];\)\)

Note, that, in out 1D model, pure heat pumping operation means not only the
            isothermal outside condition T2 = T1, but also adiabatic boundary conditions
            over the length of the pellet  0<x<Lprobe !

Show[pl1, pl2];

Absorbed thermal power

To match a practically relevant situation, a constant cross-sectional area  Aq = π/4*d^2 = 10 mm^2
            is assumed  in this section ;  Aq  is added to the set of parameters "param2"  :

(*assume cross - 
          sectional area :   Aq  =  10 mm^2  = 
          0.00001 m^2     *)  

       param2 = 
      Append[param, Aq -> 0.00001]

The absorbed thermal power  P=jq(0)*Aq  will be discussed now as a function of the temperature difference dT
            and the electric current  IA = jel * Aq  :

P[dT_, IA_] = 
      Aq*jqx0   /.  jel -> IA/Aq    // Simplify

Discussion of parameter's influence :  3DPlot of  the absorbed thermal power  P(dT, IA) :

\(\(\(\ \)\(Plot3D[
        P[dT, IA]\  /. \ param2, \ {dT, 0.1, 70}, \ {IA, 0, 
          100}, 
	\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ PlotPoints -> 
          50, \ PlotRange -> {0, 
            1.5}, 
	\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ BoxRatios -> {2, 
            2, 1}, \ ViewPoint -> {\(-2\), \(-2\), 
            1}, 
	\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ AxesLabel -> {" dT", "I   ", " P[dT,I] "}\ ]\ ;\)\)\)

Obviously, there is an optimum current that maximizes the Peltier effect.

            In addition to the Plot3D :  2D-plots with variations of parameters :
                      first plot :  P(dT) with electric current as parameter

\(\(\(\ \)\(\(\(Ppar[dT_, 
          IA_] = \(Aq*jqx0\ \  /. \ jel -> IA/Aq\)\ \ \  /. \ param2\  // 
          Simplify;\)\(
\)
    \)\(
\)
    \(\ \)\(Plot[\ {\ Ppar[dT, 4], \ Ppar[dT, 8], \ Ppar[dT, 12], 
          Ppar[dT, 20], 
		\ \ \ \ \ \ \ \ \ \ \ \ Ppar[dT, 30], 
          Ppar[dT, 40]}, 
	\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {dT, 
          0, 70}, PlotRange -> \ {0, 
            1.6}, 
			\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ AxesLabel -> {"  dT", "    P[dT]"}\ ]\ ;\)\(
\)
    \(	\ \ \ \ \ \ \ \ \ \ \ \ \)\)\)\)

second plot : P(IA) for various temperature differences

\(\(Plot[{\ Ppar[\ \ 0, IA], 
          Ppar[10, IA], 
		\ \ \ \ \ \ \ \ \ \ \ \ \ \ Ppar[20, IA], 
          Ppar[30, IA]\ , 
          Ppar[40, IA]\ }, 
	\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {IA, 0, 100}, 
        PlotRange -> {0, 
            1.5}, \ 
		\ \ \ \ \ \ \ \ \ \ \ \ \ AxesLabel -> {"  I", "    P[I]"}];\)\)

The curve's maximum changes linear in the last graphics;
           the optimum current is largest in the case of isothermal heat pumping at T1=T2 !

           Thus, within a broad field of operating conditions, there are two very particular cases :
               ◼   maximum cooling dT -> dTmax  at zero heat absorption ,
               ◼   pure heat pumping at isothermal outside condition.

           Values of optimum electric current and current density for these two cases

jelopti   /.  param
IAopti  = jelopti*Aq  /.  param2
jelopti2   /.  param
IAopti2  = jelopti2*Aq /.  param2

Summary: The optimum current  increases linear 
              from  I =  32.65  A  ( case  maximum cooling  jq=0  :  dTmax = 60.76 K  ) 
                 to   I =  40.95  A  ( case  isothermal pumping   jq->Max = Pmax/Aq  at  T1=T2  and dT=0, respectively).

Calculate maximum value of absorbed thermal power at isothermal outside condition T1=T2  :

Pmax = P[0, IAopti2 ]    /.  param2

Finally:  Contour-Plot of cooling power  P[dT,IA] in the range P > 0 and  dT > 0 :

\(\(\(\ \)\(\(\(MaxCurve = Round[10*Pmax]\)\(
\)
    \)

    \ \(ContourPlot[
        P[dT, IA]\  /. \ param2, {IA, 1, 100}, \ {dT, 0.1, 
          70}, \ 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 	Contours -> MaxCurve, 
        ContourShading -> 
          False, 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ ContourSmoothing -> 
          Automatic, 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ PlotRange -> {0, Pmax}, 
        AxesLabel -> {"jel", "jq"}\ ];\)\)\)\)

Isoline differences are  100  mW  cooling power  (see parameter MaxCurve !)

Maximum Peltier cooling : Parameter relations

This case has been discussed in detail in the last section :

                                               jel = ((σ  S)/L) T1 ,       T1min < T1 < Th
                                                                                   
            i.e. the optimum current density increases linear from maximum cooling to
            isothermal outside conditions :

σ*S/L*Th   /.  param    // N
σ*S/L*T1min   /.  param    // N
\(\(\(Plot[σ*S/L*Tx\  /. \ param\ , {Tx, Tc, Th}, 
        PlotRange -> {0, 5000000}, 
        AxesLabel -> {"T", \*""(Subscript[j, el])""}];\)\(
\)
    \)\)

Furthermore, in the case of maximum Peltier cooling, there is another common relation,
             which can easily be proved :
                                                                        dTmax =(σ/(2 κ))(S^2)(T1min^2),
                                                         
             i.e. the maximum temperature difference is always proportional to the square of the
             cold side temperature.

dTmax

    
σ*S^2/2/κ*T1min^2  // 
      FullSimplify

    
% - %%   /. Z -> (((S^2) σ)/κ)  // 
      FullSimplify

Coefficient of performance (C.O.P.)

Definition:   coefficient of performance    η = ((jq(x = 0))/(jel *  U))     with    U = | φel (x=L) | ,
           
           i.e. the coefficient of performance (C.O.P.) is defined as quotient of  pumped heat 
           from cold side x=0  through Aq per second to the electric power supplied 
           to the system of cross-section Aq

Define C.O.P. and plot as a function of cold side temperature and electric current density

Clear[T1, η];

    

    φel[L]   /. param 

    

    jq0[jel, T1]  /.  param

η[jel_, T1_] = 
        jq0[jel, T1]/(jel*Abs[φel[L]]) /. param ;

3D-plot  (note that  η > 1  is possible for heat pumps !)

\(\(\(\ \)\(Plot3D[η[jel, T1], {jel, 1, 2*10\^6}, {T1, 300, 
          250}, 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ ViewPoint -> {1, \ \(-2.5\), 1}, 
        PlotPoints -> 
          40, 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ AxesLabel -> {\*""(Subscript[j, el])"", \*""(Subscript[T, 1])"", "η"}];\)\)\)

2D-plots for a given cold side temperature

\(\(Plot[η[jel, 295], {jel, 100000, 5*10\^6}, 
        AxesLabel -> {\*""(Subscript[j, el])"", "η"}, 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ PlotRange -> {{0, 
              5*10\^6}, {\(-1\), 8}}, 
        AxesOrigin -> {0, 0}];\)\)

\(\(Plot[η[jel, 290], {jel, 100000, 5*10\^6}, 
        AxesLabel -> {\*""  (Subscript[j, el])"", "η"}, 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ PlotRange -> {{0, 
              5*10\^6}, {\(-1\), 5}}, 
        AxesOrigin -> {0, 0}];\)\)

\(\(Plot[η[jel, 280], {jel, 100000, 5*10\^6}, 
        AxesLabel -> {\*""  (Subscript[j, el])"", "η"}, 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ PlotRange -> {{0, 
              5*10\^6}, {\(-1\), 2}}, 
        AxesOrigin -> {0, 0}];\)\)

Important ! 
           The current for maximum C.O.P.  is different from the optimum current for maximum heat transport !
          
           We show this for the last graphics (case T1=280 K)  
         (for other temperatures: please change parameters "maxeta" and "maxheat"  if necessary !)

Tc = 280;

    
maxeta = 2; 
maxheat = 120000;

    
\(\(Show[
        GraphicsArray[
\ \ \ \ \ \ {Plot[η[jel, 
                Tc], {jel, 100000, 5*10\^6}, 
              AxesLabel -> {\*""(Subscript[j, el])"", "η"}, 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ PlotRange -> {{0, 5*10\^6}, {\(-0.5\), maxeta}}, 
              AxesOrigin -> {0, 
                  0}, 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ DisplayFunction -> 
                Identity], 
\ \ \ \ \ \ Plot[
              jq0[jel, Tc], {jel, 15000, 5*10\^6}, 
              AxesLabel -> {\*""(Subscript[j, el])"", \*""(Subscript[j, q])""}, 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ PlotRange -> {{0, 5*10\^6}, {\(-20000\), maxheat}}, 
              AxesOrigin -> {0, 
                  0}, \ \ \ \ \ \ \ \ \ \ \ 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ DisplayFunction -> Identity]}], 
        DisplayFunction -> $DisplayFunction, 
        ImageSize -> {600, 
            170}];\)

 (*maximum\ of\ C . 
            O . P . \((jel)\) - 
        curve\ \ \ *) 
\)
\(\ \ \ \ \ \ \ \ \ \ \ \(\(Print["maximum of C.O.P.(jel)-curve at jel : calculate e.g. with bisection method !"\ \ ]\)\(
\)\(
\)
    \( (*maximum\ of\ C . O . P . \((jel)\) - 
        curve\ \ \ *) \)\(\ \ \ \ \ \ \ \ \)\)\)
           Print["maximum of pumped heat(jel)-curve at \
jel = ",      σ *S/L*Tc /. 
        param    ]

Data-Table of C.O.P.(jel)  is the simplest way  to approximatly determine to maximum of  the  C.O.P.(jel)-curve

TableForm[
      Table[{jel, η[jel, Tc]}, {jel, 500000, 1200000, 
          50000}]]

Special case : Maximum heat transport
       Calculation of C.O.P.  at optimum current  for a given cold side temperature
       ( all outputs including T(x),  jq(x)  and  phiel(x)  along the pellet )

Tc  = 
      280                         

    
dT = T2  - Tc    /.  param

      
jeloptTc  = σ *S/L*Tc /.  param

      
T[x_] = -((((jeloptTc))^2)/(2  σ κ)) (x^2) + 
            x (1/L) (dT + ((((jeloptTc))^2)/(2 σ κ)) (L^2)) + Tc   /. param   // N
Plot[T[x]  /.  param, {x, 0, Lprobe}, 
        AxesLabel ->  {"x", "T(x)"}, 
        PlotRange -> {200, 320}];

              
 T[0]
jq[x_]  = jeloptTc*S*T2 - 
          dT ((κ/L) + 
                S*jeloptTc) - ((((jeloptTc))^2)/(2 σ)) 
            L /.param 
Plot[jq[x]  /.  param, {x, 0, Lprobe}, 
        AxesLabel ->  {"x", "jq(x)"}, 
        PlotRange -> {0, 140000}];

    
H = jq[0]
phiel[x_] = 
        S ((((jeloptTc))^2)/(2  σ κ)) (x^2) - 
            x ((jeloptTc/σ) + 
                  S (1/L) (dT + ((((jeloptTc))^2)/(2 σ κ)) (L^2)))   /. param  ;
\(\(Plot[phiel[x], {x, 0, Lprobe}, 
        AxesLabel -> {"x", \*""\!\(φ\_el\)(x)""}, 
        PlotRange -> {\(-0.2\), 0}];\)\)
U = Abs[phiel[L]]   /.  param

      
COP =  H/jeloptTc/
        U

  

    
η[jel, 
        Tc]   /.  jel->jeloptTc

Internal resistance

Under isothermal conditions, the internal resistence of the active material is  (Subscript[R, i]) = L/(σ*Aq) ; 
           thus, the voltage drop (electric potential difference) is  U=I*(Subscript[R, i]) =  L / σ * jel 
           We compare this result with the theoretically calculated potential difference  Δφel = φel(x=L). 
           Both functions are plotted as functions of the electric current density jel :

Tc  = 
      280                     

    
dT = T2  - Tc    /.  param

      
 Uinternal[jel_] = L/σ*jel    /.  
        param

      
phielL[
        jel_] = S ((((jel))^2)/(2  σ κ)) (x^2) - 
            x ((jel/σ) + 
                  S (1/L) (dT + ((((jel))^2)/(2 σ κ)) (L^2)))   /. x -> L   /. 
        param 

      
\(\(Plot[{Abs[phielL[jel]], Uinternal[jel]}, {jel, 0, 
          3000000}, 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ PlotStyle -> {RGBColor[1, 0, 1], RGBColor[0, 0, 1]}, 
        AxesLabel -> {\*"" (Subscript[j, el])"", \*""(Subscript[φ, el])(x=L)""}];\)\)

red line :  voltage drop (electric potential difference) U=I*(Subscript[R, i])  ;  
                              blue line :  theoretically calculated potential difference  φel(x=L)

It can easily be proved that both functions are identical in the case of isothermal outside conditions T1=T2  !

            Consequence for the coefficient of performance    η = ((jq(x = 0))/(jel *  U))  <   (Subscript[η, isotherm])= (σ/L) ((jq(x = 0))/(jel^2))    !
                
             Example: Tc=290 K ( change  Tc  and test  the  C.O.P. curve  ! )

\(\(\(\ \)\(\(\(Clear[Tc];\)\(
\)
    \)

    \ \(\(Tc\  = 
      290\)\(\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \)\( (*\ \ \ note\  : \ \ C . O . P . \  > \ 0\ \ only\ if\ \ Tc\  > 
          T1min\ \ \ *) \)\(
\)
    \)

    dT = \(\(T2\  - Tc\ \ \  /. \ param\)\(
\)
      \)

    \ η[jel, T1]\  /. \ \(\(T1 -> Tc\)\(
\)
      \)

    \ \(\(jq0[jel, T2]\  /. param\)\(
\)
    \)

    \ \(\(ηisotherm[jel_] = 
      jq0[jel, T2]\ /jel^2*σ/L\ \  /. param\)\(
\)
    \)

    \ \(Plot[{ηisotherm[jel], η[jel, Tc]}, {jel, 0, 
          5000000}, 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ PlotStyle -> {RGBColor[1, 0, 1], 
            RGBColor[0, 0, 
              1]}, 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ PlotRange -> {\(-1\), 5}, 
        AxesLabel -> {\*""  (Subscript[j, el])"", "η"}];\)\)\)\)

red line :  isothermal C.O.P.  ;   blue line : theoretically calculated  C.O.P

... to be continued ......