One-dimensional  Modelling of
                 Thermoelectric  Cooling  (II)

Dr. W. Seifert,  S. Walczak,  
           FB Physik, Department of Theoretical Physics, Martin-Luther-University Halle-Wittenberg, Germany

Introduction

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)
              Concerning the theoretical description of thermoelectric field coupling we refer
              to the previous notebook "PeltierCooling.nb"  and to our paper published in
              physica status solidi(a) 194, 1, 277-290 (2002).
             
              This notebook presents an introduction to thermoelectric effects (Seebeck, 
              Peltier, Thomson; see e.g.: www.thermoelectrics.com/introduction.html)
              within the frame of an one-dimensional approach assuming constant material
              properties and steady-state conditions. 
              Besides, the influence of Thomson's heat on the heat equation will be discussed.
              
              Assuming 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 the 
              principles of conservation of charge and energy        
					∇· (Overscript[jel, →]) = 0
					∇· (Overscript[jq, →]) =  (Overscript[jel, →]) ·(Overscript[E, →])					(2)
              and  the constitutive relations (1). 
              They are :
			 (∂/∂)x(-κ ∂(T/∂)x+ S T jel) = ((jel^2)/σ) + S jel  ∂(T/∂)x			(3)
				E ≡ -  ((∂(Subscript[φ, el])(x))/(∂x)) =  (jel/σ)+ S ∂(T/∂)x				(4)

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

General Hint:
	We use SI-units for all parameters, e.g:
		- the length L of the TE-active material in m,
		- the temperature T in K,
		- the Seebeck-coefficient S in V/K,
		- the thermal conductivity κ in W/(K m)
		- the electrical conductivity σ in (Ω ((m))^(-1))
		- the electric current jel in A/(m^2)
		- the thermal flux jq in W/(m^2)

Input  parameters  (SI-units) :  
                                        Lprobe :  length of the active material (typical pellet length Lprobe = 1...5 mm)
                                                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};
 param = 
      Append[param, 
        Z -> (((S^2) σ)/κ) /.  param]

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 !

Seebeck effect

The Seebeck effect  (Thomas Johann Seebeck,  
            1821)  is converting 
            heat differences directy to electricity.

      We specialize  Eqn .   (1)   for the case      (Subscript[j, el])=0    and     ∇T≠0   :

                    (Overscript[jel, →]) =σ(Overscript[E, →])-σ S ∇T          ⟹           ∇(Subscript[φ, el])+S ∇T= 0.

             From     (Overscript[E, ⇀])=-∇(Subscript[φ, el, ])  we  get               U=\(-\(∫\_\(P\_1\)\^P\_2∇φ\_el\ ·dr\&⇀\)\)      
                 

             \("Supposing S≡const., we finally find  for the "\ "thermoelectricvoltage"\)

                                 U=\(∫\_\(P\_1\)\^P\_2 
                S\ \ ∇T\ ·\ dr\&⇀\) = S (Subscript[((∫_((Subscript[P, 1])))^P), 2])dT =  S  ΔT.   

The Seebeck  coefficient   S  is defined as the open circuit voltage 
\(\(\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \)\(produced\ \ between\ two\ point\ on\ a\ conductor, \ where\ an\ uniform\)\(\ \)\)
\(\(\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \)\(temperature\ \ difference\ of\ \ \ 1\ K\ \ exits\ between\ those\ \ points . \ 

\ \ \ \ \ \ \ \ \ \ \ The\ Seebeck\ effect\ is\ the\ reverse\ of\ the\ Peltier\ effect\ and\ \ the\ basis\ of\ thermoelectric\)\)
             power generation.  ( see  e . g .   Phys .  Rev .  Lett .  91,  
              148301 (2003) ).
        This problem will be treated in a seperate notebook.

Peltier effect

The peltier effect occurs when you take any two members of the thermoelectric 
           series and connect wires made of them to form a circuit with two junctions.
           In the presense of a current a temperature difference between the juctions is produced.

            The Peltier effect (discovered by the physicist  Jean Charles Athanase Peltier in 1834)
            is the reverse of the Seebeck effect and the basis of thermoelectric heat
            transport / thermoelectric cooling.

1D Constant Properties Model (1D CPM)

Assuming the simplifying case of constant material parameters  S, κ, σ, 
            Eq. (3) reduces to the simple differential equation 
           
                              (T^″)(x) = - (Subscript[c, o])  ;     (Subscript[c, o])=  ((jel^2)/(σ κ))				(5)

This differential equation has been investigated in detail in the previous  
            notebook "pss194-1.nb". Here, we recall only some results that are 
            necessary for a comparison of  Peltier effect and Thomson effect.
            
            Hint: We only consider boundary conditions of  first kind  ( T(0)=T1, T(L)=T2 ) .
                      Mixed boundary conditions  T(L)=T2, and heat flux  jq(0)=0  
                      are also possible, but will not be considered further here.

deq1 = T''[
          x]==-c0

deqsol = DSolve[{deq1, T[0] == T1, T[L] == T2}, T[x], x] // 
          Simplify // Flatten

notation:  TPe(x)  for  Peltier effect temperature profile

TPe[x_] = T[x] /. deqsol[[1]]

test boundary conditions:

TPe[0]
TPe[L]

calculate electric potential  φel(x)  
            notation:  phiPe  for  Peltier effect  electric potential

diffeqphiPe = -D[phi[x], x] == 
        (jel/σ) + S*D[TPe[x], x]

solphiPe = 
      DSolve[{diffeqphiPe, phi[0] == 0}, phi[x], x] // Flatten

phiPe[x_] = phi[x] /. solphiPe [[1]]

prove   φel(0)=0

phiPe[0]

Calculate heat flux  jqPe as a function of x

jqPe = -κ D[TPe[x], x] + S*jel*TPe[x]

Calculate total heat flux at the cold side   (Subscript[j, q])(x = 0)

jqPe /. x -> 0

Define this heat flux as function of the temperature difference dT = T2-T1  and the electric current density  jel :

jqPe0[dT_, jel_]= jqPe /. x -> 0  /.  T1 -> T2 - dT   /. c0  -> jel^2/κ /σ  //Simplify

Example: for a defined value of  jel = 1000000 :  plot  jq(dT) 

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

\(\(plot1 = 
        Plot[jqPe0[dT, 1000000] /. param, {dT, \ 0, 
            50}, \ \ AxesLabel -> \ {"ΔT", \*""(Subscript[j, q])(ΔT)""}];\)\)

3D plot of  jq(dT,jel)⩾ 0

\(\(Plot3D[
        jqPe0[dT, jel] /. param, {jel, 0, 9*10\^6}, {dT, 0, 
          70}, 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ PlotRange -> {0, 150000}, PlotPoints -> 25, 
        ViewPoint -> {1.7, \ \(-2.1\), \ 
            0.9}, 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ AxesLabel -> {\*""(Subscript[j, el])"", "ΔT", \*""(Subscript[j, q])(ΔT,(Subscript[j, el]))   ""}, 
        ImageSize -> {399, 222}];\)\)

Within a broad field of operating conditions, there are two very particular cases :
                      •  the curved line  jq  = 0  represents the case of maximum cooling  (T1 = T1,min ),
                      •  the curved line  dT =0  represents the case of maximum heat pumping 
                         at isothermal outside conditions  T1 = T2 !
            Obviously, there is an optimum current that maximizes the Peltier effect.        
            (For parameter variations and a systematic discussion of operation conditions : see "pss194-1.nb"  )

Example: 2D-plots with variations of parameter :  jq (jel) with dT as parameter

\(\(fig1 = \ 
        Plot[{jqPe0[0, jel] /. param, jqPe0[10, jel] /. param, 
            jqPe0[20, jel] /. 
              param, \ \ \ \ 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ jqPe0[30, jel] /. param, 
            jqPe0[40, jel] /. 
              param}, 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {jel, 0, 10*10\^6}, 
          PlotRange -> {0, 
              150000}, 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ PlotStyle -> RGBColor[1, 0, 1], 
          AxesLabel -> {\*""(Subscript[j, el])"", \*""(Subscript[j, q])""}];\)\)

Isothermal heat pumping

In this section, isothermal heat pumping ( i.e.  (Subscript[T, 1]) =  (Subscript[T, 2])) will be investigated.

Tiso[x_] = 
      TPe[x] /. {T1 ->  T2, 
            c0 ->  ((jel^2)/(κ σ))} // 
        FullSimplify

The related heat flux is

jqiso[x_] = 
      jqPe /. {T1 ->  T2, c0 ->  ((jel^2)/(κ σ))} // 
        FullSimplify

Calculate optimum electric current density

eq1tpy = 
      Solve[D[jqiso[0], jel] == 0, jel] // 
        Flatten
jelisoopti = jel /. eq1tpy[[1]];
jelisoopti  /. param

Calculate and plot temperature profile at optimum electric current density

Tisopti[x_] = 
        Tiso[x] /. jel -> jelisoopti;
\(\(bild1 = 
        Plot[Tisopti[x] /. param, {x, 0, 
            Lprobe}, 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ AxesLabel -> {"x", "T(x)"}, 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ PlotStyle -> {{RGBColor[1, 0, 
                  1]}}];\)\)

Note, that, in our 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 !

Example: calculate and plot  electric potential under isothermal outside conditions

phiiso[x_] = 
      phiPe[x] //. {T1 -> T2, 
            c0 ->  ((jel^2)/(κ σ)), 
            jel -> jelisoopti} // Simplify
\(\ \(bildpot1 = 
        Plot[phiiso[x] /. param, {x, 0, 
            Lprobe}, 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ AxesLabel -> {"x", \*""(Subscript[φ, el])""}, 
          PlotStyle -> {RGBColor[1, 0, 1]}];\)\)

Maximum electric potential difference (voltage drop)

phiiso[Lprobe] /. param // N

Maximum Peltier cooling

Arbitrary combinations of parameters are not usefull to understand thermoelectric cooling.
            A systematic investigation of maximum Peltier Cooling is presented in the notebook "pss194-1.nb". 
            As to the case of maximum cooling, we recall here only results that are important for comparison
            to the Thomson effect.

Calculate the optimum current for maximum cooling

eq2tpy = Simplify[Solve[D[jqPe0[dT, jel], jel] == 0, jel], 
          dT == 
            T2 - T1]//Flatten

jelmkopti = jel /. eq2tpy[[1]];

Note, that the equation    jopt = ((σ  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 !

Strategy: Calculate heat flux jq(x=0) at the cold side 
                             and evaluate optimum current and maximum temperature difference

eq3tpy = Solve[jqPe0[dT, jel] == 0, dT] // Simplify // 
        Flatten

deltaT = dT /. eq3tpy [[1]]

eq4tpy = FullSimplify[
          Solve[D[deltaT, jel] == 0, jel], κ > 
            0]//Flatten

eq4tpy /. param

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

jelmkopti = jel /. eq4tpy[[2]]
%   /.  param

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

dTmax = FullSimplify[
          deltaT /. jel -> 
              jelmkopti]//Expand

T1mkmin = T2 - dTmax

values ...

dTmax /. param
T1mkmin /. param

Plot temperature profile in the case of  maximum cooling

Tmk[
        x_] = TPe[x] /. {c0 ->  ((jel^2)/(κ σ)), 
              T1 -> T1mkmin} /. jel -> jelmkopti // 
        FullSimplify
\(\(bild2 = 
        Plot[Tmk[x] /. param, {x, 0, Lprobe}, 
          AxesLabel -> {"x", "T(x)"}, \ \ \ \ \ 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ PlotStyle -> RGBColor[1, 0.2, 0]]\ ;\)
\)

It is interesting to note that the temperature gradient vanishes at the warm side :  ((∂T )/∂)x(x=L)= 0 .
            This is a consequence of the boundary condition  T( x=L )=T2 = const. !

D[Tmk[x], x] /. x -> Lprobe    /.  param   // 
      Simplify

Example: calculate and plot electric potential in the case of maximum cooling

phimk[x_] = 
      phiPe[x]  //. {c0 ->  ((jel^2)/(κ σ)), 
            T1 -> T1mkmin, jel -> jelmkopti} // 
        Simplify
\(\(bildpot2 = 
        Plot[phimk[x] /. param, {x, 0, Lprobe}, 
          AxesLabel -> {"x", \*""(Subscript[φ, el])""}, 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ PlotStyle -> RGBColor[1, 0.2, 0]]\ ;\)\)

Maximum electric potential difference (voltage drop)

phimk[Lprobe]  /. param

Summary: Comparison of both cases under discussion

Comparison of temperature profiles 
                                                lila line :  isothermal outside conditions
                                   red dashed line :  case of  maximum cooling

bild3 = 
        Show[bild1, bild2, PlotRange -> {220, 330}];

Comparison of calculated elctric potential
                                             lila line :  isothermal outside conditions
                                red dashed line :  case of  maximum cooling

bildpot3 = Show[bildpot1, bildpot2];

Small differences in the electric potential profile (Subscript[φ, el])(x) are due to the
           different temperature profiles inside the active material ...

Thomson effect

William Thomson (later Lord Kelvin) established the relationship between the Peltier
          and Seebeck coefficients :  Π = T· S ,  and predicted the third thermoelectric effect,
          the Thomson effect (1854):  if a temperature difference exists between any two 
          points of a current-carrying conductor, heat is either evolved or absorbed depending 
          upon the material.

1D Quasi-CPM

Suppose that a current passes along an active thermoelectric material in the presence 
          of a temperature gradient.
          Any  theoretical decription of  the Thomson effect starts with the assumption that 
          the Seebeck coefficient S  is a  temperature dependent material property. 
          With  Π(T) = T · S(T)  we get for the gradient
         
                             ∇ Π = d(Π/(dT )) ∇T = ( S(T) + T (dS/(dT )))  ∇T 
         
          Introducing the Thomson coefficient  (Subscript[μ, E]) by one of the definitions
         
                               (Subscript[μ, E]) :=   d(Π/(dT )) - S(T)  =  T (d/(dT ))( d(Π/(dT )))  =  T (dS/(dT ))   ,                         (6)
                               
          we obtain the following relations : 
                            
                                ∇ Π =  ( (Subscript[μ, E])+ S ) ∇T           and        S (T) =  (∫_0^T)(( (Subscript[μ, E])(
                T'))/T)' dT'  .
                                
          Thus, the temperature dependence of  (Subscript[μ, E])(T) and  S(T)  is correlated ;  
          both coefficients S and  (Subscript[μ, E]) can be evaluated e.g. from experimental available  Seebeck data
          if the real S(T) dependence can be approximated by a logarithmic  S(T) law.
          
          Within the frame of a quasi-CPM, we assume now  S(T), but  (Subscript[μ, E]), κ, σ = const. , 
          so that the real S(T) dependence is approximated by a logarithmic S(T) law :
          
                                        S(T) = (Subscript[S, ref]) + (Subscript[μ, E])ln (T / (Subscript[T, ref]))  .
                                        
           The (1D) differential equation to be solved is 
           
                                     T''(x)  - (Subscript[k, o]) T'(x)  = -(Subscript[c, o])    with   (Subscript[k, o]) = (Subscript[μ, E])\(TraditionalForm\`\^j\_el\%\(\(\ \ \ \)\(o\)\)\/κ\)
                                       
            Investigating the Thomson effect, several points of view are possible :
          
             (i)  Compare the deq's solution T(x)  both for  (Subscript[μ, E])= 0  and  (Subscript[μ, E]) ≠ 0   at the same BC's 
                   and the same electric current density jel
                    Problem:  we do not consider different optimum conditions  for  (Subscript[μ, E]) = 0  and  (Subscript[μ, E])≠ 0 
                                    which lead to different optimum jel values ...
             
            (ii)  Compare CPM and Quasi-CPM 
                   Problem:  for a comparison between the more realistic case considering the Thomson effect 
                                    and the simplifying case of constant material properties (CPM), the central question 
                                    is the appropriate choise of the free parameters  (Subscript[S, ref])  and   (Subscript[T, ref]) .
                    Detailed investigations have shown that best accordance has been achieved assuming  S 
                    (within the CPM)  to be the volume average <S>,  i.e. in our 1D model :  
                                                                              S =  <S> =   1 / L (∫_0^L) S( 
            T(x) ) dx      
                    Note that the volume average is influenced  by the temperature profile T(x) ,
                    i.e. a self-consistent procedure is necessary to find the appropriate  <S>  !
                    
            In this notebook, we refer primarily to case (i). Detailed investigations about case (ii) 
            was presented in :     physica status solidi(a) 194, 1, 277-290 (2002).
     
            The following strategy is close to that used in our previous notebook "pss194-1.nb"

Quasi-CPM and boundary conditions of first kind

We assume now that the real  S(T) dependence can be approximated by a 
            logarithmic S(T) law together with a constant value of  (Subscript[μ, e])= 120 μV/K
            which can be considered as a representative, referenced value for
            non-degenerate extrinsic semiconductors.
            
            Besides, we choose for further calculation  :      (Subscript[S, ref]) = S      and       (Subscript[T, ref]) = Th  .
                                
            We compare now  S(T) with the given constant  Seebeck coefficient  S = 210*(10^(-6)) V / K  :

param2 = {param, {μ -> 0.12*(10^(-3)), 
            Tref -> 
              Th}}   //Flatten

   SvonT[
        T_] = Sref +  μ*Log[T/Tref ]    /.  
          Sref ->  
            S    /. param2

\(\ \ \ \(plotSvonT = 
        Plot[{S\  /. \ param2, SvonT[T]}, {T, 10, 400}, 
          AxesLabel -> {"T", "S(T)"}, 
          PlotRange -> {0, 
              0.0003}, 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ PlotStyle -> {RGBColor[1, 0, 1], 
              RGBColor[0.6, 0.2, 0.5]}]\ ;\)\)

Proof :  calculate Thomson coefficient

T*D[SvonT[T], T]

Note  that  SvonT(T) ⩽  S  for  T ⩽ T2 = 300 K ;  from this we conclude that  

                            <S> =   1 / L (∫_0^L) S( T(x) )dx    will be certainly smaller than S !   
                                        
             (   (Subscript[S, ref])≃ 223*(10^(-6))  V/K   would  lead  to    <S(T)> = S = 210*(10^(-6))   V/K   )
                     
            Thus, an  immediate comparison  of  CPM and  Quasi-CPM  is pretty difficult ...
            To check this, the differential equation including Thomson's heat has to be solved :
           
                 T``(x) -  (Subscript[k, o]) T`(x) =  - (Subscript[c, o])  with   (Subscript[k, o])=(Subscript[μ, e])(((Subscript[j, el])^o)/κ)  and   (Subscript[c, o])= ((((((Subscript[j, el])^o)))^2)/(σ κ))   .          (7)
                     
            Eq. (7) can be evaluated from Eq. (3). 
            Note, that, for  (Subscript[μ, E])= 0  and  (Subscript[k, o])= 0,  Eq.(5) holds.

deq2 = 
      DSolve[{T''[x] - k0 T'[x] == -c0, T[0] == T1, 
            T[L] == T2}, T[x], x] // Flatten

notation:  Tho(x)  for  Thomson effect  temperature profile

Tho[x_] = T[x] /. deq2  [[1]]

Note:  The temperature profile depends on (Subscript[μ, E]) ;  
                        but the absorbed thermal flux also depends on  S(T)  ...

Define the spatial  Seebeck coefficient  S(x) = S(T(x)) = S + (Subscript[μ, e])Log(((T(x))/(Subscript[T, ref])))

varS[x_] = 
      S + μ*Log[Tho[x]/Tref]    /. Tref -> T2      // 
        Simplify

Further strategy now as used  in the notebook "pss194-1.nb" :
            Define the thermal flux  as a function of jel  and the temperature difference dT=T2-T1 :

jqTho[dT_, 
        jel_] = -κ D[Tho[x], x] + 
            Tho[x]*varS[x]*
              jel 
                                         //. {c0 -> 
              ((jel^2)/(σ κ)), 
            k0 -> ((μ*jel)/κ), T1 -> T2 - dT}   // 
        FullSimplify

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

jqTho0[dT_, jel_] = 
      jqTho[dT, jel] /. x -> 0  // FullSimplify

Plot heat flux as function of the temperature difference dT=T2-T1 for an 
            electric current density  jel =1 000 000 (A/(m^2))

\(\(plot2 = 
        Plot[jqTho0[dT, 1000000] /. param2, {dT, 0, 
            50}, 
\ \ \ \ \ \ \ \ \ \ \ \ \ PlotStyle -> RGBColor[1, 0.4, 0], 
          AxesLabel -> {"ΔT", \*""(Subscript[j, q])(ΔT)""}];\)\)

Comparison with Peltier effect  ( black solid line : (Subscript[μ, E ])=0 )

Show[plot1, plot2, 
        PlotRange -> {-10000, 60000}];

Two very particular cases - already known from the Peltier effect - arise here modified by Thomson's heat :
                         •   maximum cooling dT -> dTmax  at zero heat absorption ,
                         •   pure heat pumping at isothermal outside condition.
                
            3DPlot of  the absorbed thermal power  jq(dT, jel) :

\(\(Plot3D[
        jqTho0[dT, jel] /. param2, {jel, 10\^\(-3\), 9*10\^6}, {dT, 0, 
          70}, 
PlotRange -> {0, 150000}, PlotPoints -> 25, 
        ViewPoint -> {1.7, \ \(-2.1\), \ 
            0.9}, 
\ \ \ \ AxesLabel -> {\*""(Subscript[j, el])"", "ΔT", \*""(Subscript[j, q])(ΔT,(Subscript[j, el]))   ""}, 
        ImageSize -> {399, 222}];\)\)

The figure is analoguous to the Peltier effect  :
                      •  the curved line   jq = 0  represents the case of maximum cooling  (T1 = T1min ),
                      •  the curved line  dT = 0  represents the case of maximum heat pumping 
                         at isothermal outside conditions  T1 = T2 !

Parameter variations : 2D-plot of  jq(x=0, dT)  with electric current as parameter

\(\(Plot[{jqTho0[dT, 100000] /. param2, 
          jqTho0[dT, 500000] /. 
            param2, 
\ \ \ \ \ \ \ jqTho0[dT, 1000000] /. 
            param2, jqTho0[dT, 1500000] /. param2}, {dT, 0, 
          60}, \ \ \ \ 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ PlotRange -> {0, 100000}, 
        AxesLabel -> {"dT", \*""(Subscript[j, g])""}];\)\)

Parameter variations :   jq(x=0, jel)  for various temperature differences dT

\(\(fig2 = 
        Plot[{jqTho0[0, jel] /. param2, jqTho0[10, jel] /. param2, 
            jqTho0[20, jel] /. 
              param2, 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ jqTho0[30, jel] /. param2, jqTho0[40, jel] /. param2}, {jel, 10\^\(-3\), 
            10*10\^6}, 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ PlotStyle -> RGBColor[0, 0, 1], PlotRange -> {0, 150000}, 
          AxesLabel -> {\*""(Subscript[j, el])"", \*""(Subscript[j, q])""}];\)\)

Obviously, there is also here an optimum current that maximizes the effect.
           We compare fig1 and  fig2 :

Show[fig1, fig2];

This comparison shows directly the influence of  Thomson's heat
            ( lila lines: Peltier effect  (Subscript[μ, E])= 0 ;  see fig1 above ):

Isothermal Thomson effect

An inhomogeneous temperature distribution (plotted below)
            exits also in the case of  isothermal outside conditions T1 = T2

Tthiso[x_] = 
      Tho[x] /. {T1 -> T2, c0 -> ((jel^2)/(σ κ)), 
            k0 -> ((μ*jel)/κ)} // Simplify

Calculate the related total heat flux

jqthiso[x_] = 
      jqTho[dT, jel] /. dT -> 0   // FullSimplify

We search  for the optimum current..
             Due to the non-algebraic equation,  a straight-forward calculation is not possible ...

Solve[D[jqthiso[0], jel] == 0, 
      jel]

Alternative :  Plot the partial derivative  ∂((Subscript[j, q, iso])/∂) (Subscript[j, el])

func[jel0_] = 
        D[jqthiso[0], jel] /. param2;
\(\(plot3 = 
        Plot[func[jel], {jel, 0, 6*10\^6}, 
          AxesLabel -> {\*"" (Subscript[j, el])"", \*""∂(Subscript[j, q])/ ∂(Subscript[j, el])|x=0""}];\)\)

FindRoot can help to find the optimum current ( a starting value is needed ! )

jelthisoopti = 
      jel /. FindRoot[
            func[jel] == 0, {jel, 4*(10^6)}]  [[1]]

note:  jelisoopti (Peltier effect with (Subscript[μ, E])= 0 ) and  jelthisoopti  ( (Subscript[μ, E])= 120 μV/K )  have different values  ...

jelisoopti /. param
jelthisoopti

.... that is the reason why we get different temperature profiles ...

Tthiso[
        x]                                       /. jel -> jelisoopti
Tthisoopti[x_] = 
      Tthiso[x] /. jel->jelthisoopti

Summary:  Plot of the resulting temperature profiles

\(\(bildth1 = 
        Plot[Tthisoopti[x] /. param2, {x, 0, Lprobe}, 
          AxesLabel -> {"x", "T(x)"}, 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ PlotStyle -> 
            RGBColor[0, 0, 1], 
          DisplayFunction -> 
            Identity];\)
\)
\(\(bildth3 = 
        Plot[\(Tthiso[x] /. jel -> jelisoopti\)\  /. param2, {x, 0, 
            Lprobe}, 
          AxesLabel -> {"x", "T(x)"}, \ \ \ 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ PlotStyle -> {{Dashing[{0.035}], RGBColor[0, 0, 1]}}, 
          DisplayFunction -> 
            Identity];\)
\)
Show[bild1, bildth1, bildth3, 
        DisplayFunction -> $DisplayFunction] ;

Comparison of calculated temperature profiles at isothermal outside conditions T1=T2
                                lila solid  line :     Tiso(x)          /.  jel -> jelisoopti;   i.e   (Subscript[μ, E]) = 0 ;               jel = 4.095×(10^6)      ( bild1 )
                   blue long dashed line :      Tthiso(x)       /.  jel -> jelisoopti;          (Subscript[μ, E]) = 120 μV/K ; jel = 4.095×(10^6)      ( bildth3 )
                               blue solid line :    Tthisoopti(x)  /.  jel -> jelthisoopti;       (Subscript[μ, E]) = 120 μV/K ; jel = 4.553721106364085`*^6  ( bildth1 )

Thomson effect : Maximum cooling

We consider now maximum cooling ( case  T1 < T2 ).
           The optimum curent density for maximum temperature difference   Δ(Subscript[T, max]) will be of particular interest.
           As in the previous section, a straight-forward calculation is not possible....

jqTho0[dT, jel]

Alternative: 
            Expand the logarithm into a power series !
            Here, we will only consider the first order approximation !
            
            Test: logarithmic S(T) law

SvonT[T]
Logapprox = 
      Series[SvonT[T], {T, Tref, 1}]    /.  param2   // 
        Normal

plotapprox1 = 
        Plot[Logapprox  /.  param2, {T, 0, 400}, 
          AxesLabel -> {"T", "S(T)"},  
          PlotStyle -> {Dashing[{0.015}]}];

plotapprox2 = Show[plotapprox1, plotSvonT];

Even a first order approximation (dashed line) seems to be sufficient for 
            small temperature differences around Th = 300 K  ...

Strategy now: 
            Rewrite the thermal flux jqTho0 at the cold side x = 0  (defined above)  as a function of jel  
            and the cold side temperature T1:

jqth1[T1_, jel_] = 
      jqTho0[dT, jel]    /. dT -> T2 - T1

Parameter variations:  2D-plot of  jqth1(x=0, T1) with electric current as parameter

\(\(Plot[{jqth1[T1, 500000] /. param2, 
          jqth1[T1, 600000] /. 
            param2, 
\ \ \ \ \ \ \ jqth1[T1, 1000000] /. 
            param2, jqth1[T1, 1500000] /. 
            param2}, 
\ \ \ \ \ \ \ {T1, 240, 300}, 
        PlotRange -> {\(-50000\), 100000}, 
        AxesLabel -> {\*""(Subscript[T, 1])"", \*""(Subscript[j, q])((Subscript[T, 1]))""}];\)\)

Expand  jqth1  into a power series 
             and use the first order approximation of heat flux

jtaylor = Normal[
          Series[jqth1[T1, jel], {T1, T2, 
              1}]] //Simplify

now, solving for T1 is possible ...

eqtpy6 = 
      FullSimplify[
          Solve[jtaylor == 0, T1], {μ > 0, σ > 0, T2 > 0, 
            S > 0}]  // Flatten

notation:  T1thmin1 is the minimum temperature (maximum cooling inluding Thomson's heat)

T1thmin1 = T1  /. eqtpy6

Insert T1thmin1  into  jqth1 and search  for  optimum current density ...

eqtpy7 = D[jqth1[T1, jel], jel]    /. 
            T1 -> T1thmin1    /. 
          param2;
plot4 = 
        Plot[eqtpy7, {jel, 1, 6000000}, 
          AxesLabel -> {"jel",  "jq"}, 
          PlotStyle -> RGBColor[1, 0.3, 0]];

Show[plot3, plot4];

FindRoot can help to find the optimum current

eqtpy8 = FindRoot[eqtpy7 == 0 , {jel, 3*(10^6)}]

Summary :  first order approximation results for the optimum current density  
                                to reach minimum temperature  T1thmin1

jelthoptimum1 = jel /. eqtpy8 [[1]]
T1thmin1   /. jel -> jelthoptimum1   /. param2

Note that we have different optimum current values for isothermal heat pumping and maximum cooling !

Calculate temperature profile for the case of maximum cooling

Tthmkoptimum1[
        x_] = Tho[x] //. {c0 -> ((jel^2)/(σ κ)), 
            k0 -> ((μ*jel)/κ), T1 -> T1thmin1, 
            jel -> jelthoptimum1}    /. param2

\(\ \(bildth2 = 
        Plot[Tthmkoptimum1[x] /. param2, {x, 0, Lprobe}, 
          AxesLabel -> {"x", "T(x)"}, 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ PlotRange -> {200, 300}, \ 
          PlotStyle -> {{Dashing[{0.015}], 
                RGBColor[0, 0, 
                  1]}}];\)

\  (*\ \ prepare\ plot\ for\ summary\ of\ results\ \ \ *) \ \ \)
\(\(bildth4 = 
        Plot[\(Tho[x] //. {c0 -> \((jel0)\)\^2\/\(σ\ κ\), 
                k0 -> \(μ\ jel0\)\/κ, T1 -> T1mkmin, 
                jel0 -> 
                  jelmkopti}\)
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \  /. param2, {x, 0, 
            Lprobe}, 
          AxesLabel -> {"x", "T(x)"}, 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ PlotStyle -> Hue[0.8], 
          DisplayFunction -> Identity];\)\)

Summary:  Thomson effect  ( (Subscript[μ, E]) =  120 μV/K )
                                      temperature profiles for isothermal BC's and maximum cooling
                               
                              blue solid line (bildth1):  Tthisoopti[x]  /.  jel-> jelthisoopti = 4.553721106364085`*^6;  (Subscript[μ, E])= 120 μV/K 
                 blue long dashed line  (bildth3):      Tthiso[x]     /.  jel-> jelisoopt i   = 4.095×(10^6)     ;  (Subscript[μ, E])= 120 μV/K
                   
                           blue dashed line  (bildth2):  Tthmkoptimum1[x]  /.  jel -> jelthoptimum1 = 3.3700767132565286`*^6 ;  
                                                                                                           T1 -> T1thmin1 = 243.037 K   
                       magenda solid line  (bildth4):    Tho[x]    /.  jel -> jelmkopti  = 3.26558×(10^6);   T1 -> T1mkmin = 239.236 K ,                                                                           
                                                                           note :  values   jelmkopti, T1mkmin  from  Peltier cooling  ( (Subscript[μ, E]) = 0 )

Show[bildth1, bildth2, bildth3, bildth4,  
        PlotRange -> {220, 330}, 
        DisplayFunction -> $DisplayFunction];

Finally, for the case of maximum cooling (see bild2th) :   calculate  S(x) =  S(T(x)) and the average  <S>

varS[
        x_] = S + μ*Log[Tthmkoptimum1[x]/Tref]   /.  
          param2  //Simplify

Sav = NIntegrate[varS[x], {x, 0, Lprobe}]/
        Lprobe

\(\(plotSvonx = 
        Plot[{S\  /. \ param, varS[x], Sav}, {x, 0, 
            Lprobe}, 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ PlotStyle -> {RGBColor[1, 0, 1], RGBColor[0, 0, 0], 
              RGBColor[0, 0, 
                1]}, 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ AxesOrigin -> {0, 0.00017}, 
          AxesLabel -> {"x", "  S(x)"}, 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ PlotRange -> {0.00017, 
              0.00022}];\)\)

The  average  <S>= 201.863  ( blue line )  is smaller than S ( lila line; see plotapprox2  below ! )

Comparison between Peltier effect and Thomson effect

It is also worth comparing equivalent temperature profiles with and without Thomson's heat.
            Note that the optimum current density and - in the case of thermoelectric cooling -
            the minimum temperatures may have differrent values.
             
            A)  Temperature profiles for isothermal outside conditions  (T1 = T2)
                                                   red solid line :   Peltier effect ,  (Subscript[μ, E]) = 0                        ( bild1 ) 
                                                          blue line :  Thomson effect , (Subscript[μ, E]) =  120 μV/K       ( bildth1 )

bildvi1 = Show[bild1, bildth1];

B)   Isothermal outside conditions  (T1 = T2)
                   Comparison of  the temperature profiles calculated at the same electric current density
                                              red solid line :  Peltier effect ,     (Subscript[μ, E]) = 0                        ( bild1 ) 
                                blue long dashed line :  Thomson effect , (Subscript[μ, E]) = 120 μV/K          ( bildth3 )

bildvi2 = Show[bild1, bildth3];

Due to Thomson's heat the temperature profile specifying the isothermal Thomson effect 
            (blue long dashed line) is slightly unsymmetric !

C)  Maximum cooling
                   Comparison of the temperature profiles calculated at the optimum electric current density
                                red solid line :            jel -> jelmkopti = 3.265577869445645`*^6;  T1 -> T1mkmin = 239.236  K      ( bild2 )
                          blue dashed line :    jel -> jelthoptimum1 = 3.3700767132565286`*^6 ;  T1 -> T1thmin1 = 243.037 K       ( bild2th  )

bildvmk1 = 
        Show[bild2, bildth2, PlotRange -> {230, 300}];

D)  Maximum cooling :  Influence of  Thomson's heat on T(x)
                  Comparison of the temperature profiles ( bild2, bildth4 ) calculated 
                  at the same electric current density  jel -> jelmkopti  
                  and the same minimum temperature  T1 -> T1mkmin

bildvmk2 = 
        Show[bild2, 
          bildth4, 
                            PlotRange -> {230, 300}, 
          DisplayFunction -> $DisplayFunction];