THE JOULE-THOMSON COEFFICIENT FOR A SQUARE WELL GAS

A computer algebra application in physical chemistry.

M. Hanson, 
		Department of Chemistry, 
		Humboldt State University, 
		Arcata, California 95521
		
		Bitnet: GFA001D@CCS.CSUSCC.CALSTATE.EDU
		(707) 826 3265

Introduction

There are many problems in physical chemistry that are easy to solve formally but impossibly tedious in practice.  The advent of computer algebra programs makes these problems tractable.  The determination of the Joule-Thomson coefficient μJT for a square well gas is such a problem.  One determines the enthalpy as a function of the temperature and pressure, H(T,P), differentiates with respect to P, sets the result to zero, and solves for μJT as (¶T/¶P)H.
		The three parameters of the potential (the well depth ϵ, the hard sphere size of the particle d, and the range parameter g) are set by fitting the second virial coefficient.  The resulting Joule-Thomson coefficient for argon is found to agree remarkably well with experiment.

Model

1.  A collection of the Avogadro number of particles N having the argon atom mass in volume V at an equilibrium temperature T.  

2.  The interaction between particles is given by the potential

																 		∞,    if r  ≤  d    
							U(r) =            - ϵ,    if d  <  r  ≤  g d																	(1)
															 			 0,   if g d < r
									
3.  The pressure is sufficiently low that the equation of state satisfies

												pV = RT ( 1 + B(T)/V )
																																																			(2)	
															≈  RT + B(T)p.


4.  The second virial coefficient is given by

																					 ∞								                2
										B(T)  =  2 p N  |  ( 1  -  exp(-U/kT) ) r   dr.								(3)
																			 0 
																					
For the model potential we have (see the Appendix)
																														
										B(T)  =  b  -  (a/Nϵ){exp(ϵ/kT) - 1}														(4)
										
										
	where
																	3                          3
	 				b  =  (2/3)πNd ,          a/Nϵ  =  (g   -  1)br.

Fitting the Parameters

Get and Plot the Experimental Data

ArgBofT = ReadList["argboft.list",{Number,Number}] ;
(* see NBS reference *)

gbexp = ListPlot[ArgBofT, AxesLabel->{"T/K",""},
			PlotLabel->
			"B[T]/cm^3   The Second Virial Coefficient of Argon"]

Guess d, eps(ϵ), & g, and Plot the Model Prediction

ClearAll[BofT,T,No,d,eps,g,b,gbtheor]
BofT[T_] := b*(1 - (g^3 - 1)*(Exp[eps/(k*T)] - 1) )
No = 6.0225 10^23;
d = 316.2*10^-12   (* in meters *)
eps = 69.4*k;
g = 1.85;   (* The parameters come from Guggenheim. Try others *)
b = 2*Pi*No*d^3/3 ;
gbtheor = Plot[BofT[T] 1 10^6,{T,80,700}]    (* convert to cm^3 to compare with experimental data *)

Comparison Model and Experimental Data

Show[gbexp,gbtheor]

The Energy and the Enthalpy

The energy is made up of the sum of the contributions from kinetic and potential energies  
                                            2
                 E  =   N <(1/2)ms   >  +  N <U/2>                                   (5)
                 
(Since U gives the potential for a pair of particles, U/2 gives that for a single particle.)  Using the Boltzman distribution for speeds and the Maxwell-Boltzman distribution for positions we have

                                        2
                          <(1/2)ms   >  =  (3/2)kT                                       (6)
and
                                  ∞                                 2
                <U/2>  =    |     U   N exp(-U/kT)4πr  dr.                          (7)
                              0     2   V
                              
   Using the given potential in this equation gives (see the Appendix)
   
                                                                 3
                 <U/2>  =  - (bϵ/V) exp(ϵ/kT) (g   -  1)
                                                                                                       (8)
                            =  - (a/N V) exp(ϵ/kT).
                            
    With the enthalpy given as
    
                                     H  =  E  +  pV,                                               (9)
                    
     we have
     
               H(T,p)  =  (5/2) RT + pB(T) - a p exp(ϵ/kT) / [RT + pB(T)].   (10)

The Joule-Thomson Coefficient

The Experimental Data

ArgMuJT = ReadList["argmujt.list",{Number,Number}];
plot = ListPlot[ArgMuJT, AxesLabel -> {"t/C","  "},
		PlotLabel -> 
		"          The Joule-Thomson Coefficient for Argon"
		]
		
(* see Smithsonian reference *)

The Model Prediction
(About 30 Minutes on a IIci)

ClearAll[a,R,H,k,tmp1,mujt,muJTAr,P]
a = (g^3 - 1)*b*No*eps ;
R = 8.3144 ;
H = 5 R T / 2 + P BofT[T] - a P Exp[eps/(k T)]/
																											(R T + P BofT[T])
k = R/No ;
tmp1 = Solve[D[H,P,NonConstants->{T}] == 0,
		D[T,P,NonConstants -> {T}]];
mujt[T_] = D[T, P, NonConstants -> {T} ]  /. tmp1[[1]] ;
P = 101325.(* The units are now K/Pa change them to K/atm. *)
muJTAr[T_] = mujt[T] 101325. ;

plot1 = Plot[muJTAr[t+273.15],{t,-170,300}]

Comparison Theory and Experiment

Show[plot,plot1]

(*just a spacer*)

References

Guggenheim, Edward A., Elements of the Kinetic Theory of Gases,  Oxford, 1960.

NBS Circular 564, Tables of Thermal Properties of Gases

Smithsonian Physical Tables nineth revised edition, volume 39, 1954