© 2022 Horia Metiu

WorkBookT.23.1  Calculate the coexistence curve for melting benzene

Cell 1:  Define the problem and give data.

I know the change of enthalpy and the change of volume caused by melting and the melting temperature at p=1 atm. I’ll use the Clapeyron equation
	(dp/dT) = ((Δh)/(T Δv))						(1)
to calculate how the melting temperature depends on pressure.

At p=1 atm, the melting temperature is T=278.7 K.
At this temperature and pressure the heat of melting is 
	Δh = 9.95 kJ/mol
and the change of volume upon melting is 
	Δv = 10.3 (cm^3)/mol

Since this is the only data I have I will make some approximations that are commonly made in all textbooks.  Then I can test them by comparing my calculations with the following data

	 T/deg K         p/atm
	({{311.2`,      987.`}, {393.29999999999995`,      4935.`}, {469.2`,      9870.`}})

Cell 2:  Input data

Input the heat of melting Δh, the change of volume upon melting Δv, the melting temperature Tin when p = 1 atm.

Off[Remove::rmnsm]
Remove["Global`*"] ;
Needs["Utilities`CleanSlate`"];
CleanSlate[];

Δh= 9.95 (*kJ/mol*);
Δv=10.3 (*(cm^3)/mol*);
pin=1 (*atm*);
Tin=278.7 (*K*);

The experimental values of the melting temperature at different pressures :

TExp={32.5,114.6,190.5}+278.7  ; (*K*)
pExp=0.987*{1000, 5000, 10000};(*atm *)

0.987 converted bar to atm

Display the data:

TableForm[MapThread[{#1,#2}&,{TExp,pExp}],
TableHeadings->{None,{"T K","(Subscript[p, exp]) atm"}}]

Cell 3:  Calculate the coexistence line:

• Integrating the Clapeyron equation with Δh and Δv independent of T and p gives
	p -1= ((Δh)/(Δv))ln[T/278.7] 
• The equation is telling me how the melting temperature T depends on pressure.
• Define a function p[T] that evaluates this expression. 
• Pay attention to units.  Δh is in kJ and Δv in (cm^3)/mol.
• Need Δh in liter atm and Δv in liter/mol.  
•  Use  
	1kJ = 9.8692 liter atm 
and 
	1(cm^3) = (10^(-3)) liter. 
• The calculated pressure will have units of atm.

p=.
p[T_]:=1+(9.8692* Δh/(10^-3*Δv))* Log[T/278.7](*atm*)

Calculate the pressure at the temperatures in the list TExp:

pCalc=p/@TExp  (*atm*)

Display the results and calculate the percent error (compare the calculated pressures with the measured ones):

TableForm[MapThread[{#1,#2,#3,100*(#2-#3)/#2}&,{TExp,pExp,pCalc}],
TableHeadings->{None,{"T K","(Subscript[p, exp]) atm","(Subscript[p, calc]) atm","% error"}}]

The error are very large!  The biggest error comes from neglecting the change of Δv with pressure.  The temperature does not vary a lot and the change of Δh with T should not be very large. Nevertheless, you are learning that oversimplification comes at a price.

WorkBookT.23.2  Equilibrium pressure of diamond with graphite at T = 298 K

Cell 1. Define the problem and give the data

I will use the equilibrium condition  
	μ(d; 298 K, p) = μ(g; 298 K, p)					(1)
which is the equilibrium condition.  As shown in the textbook this leads to the equation 
	p = 1- Δμ(298, 1)/Δv(298,1)						(2)
To calculate Δμ(298,1) I use chemical potential of formation
	Δμ =  (Subscript[Δμ, f])(diamond) - (Subscript[Δμ, f])(graphite)				(3)
I calculate Δv = v(d) - v(g), where v is the molar volume, from the densities.

Here is the data I need to perform the calculation of  p versus T.  From Appendix 9 I obtain
	(Subscript[Δμ, f])(diamond) = 2.9 kJ/mol
	(Subscript[Δμ, f])(graphite) = 0

I also have the densities:
	ρ(diamond)=3.51 g/(cm^3)
	ρ(graphite)=2.25 g/(cm^3)
The molar volume is 
	12 /ρ   (cm^3)/mol, 
where 12 is the molar mass of carbon.

Cell 2. Input data

Off[Remove::rmnsm]
Remove["Global`*"] ;
Needs["Utilities`CleanSlate`"];
CleanSlate[];

Δμfg=0 (*Δμ for graphite, in kJ/mol, at 298 k*);
Δμfd=2.9 (*Δμ for diamond, in kJ/mol, at 298 k*);
ρg= 2.25  (*the density of graphite under standard conditions g/(cm^3)*);
ρd=3.35 (*the density of diamond under standard conditions g/(cm^3)*);

Cell 3. Calculate Δv, Δμ and p

Δv=12*(1/ρd -1/ρg)*10^-3  (*liter/mol*)

Calculate Δμ from 
	Δμ =  (Subscript[Δμ, f])(diamond) - (Subscript[Δμ, f])(graphite);  
Use 
	1 kJ = 9.8692  liter atm

Δμ=9.8692*(Δμfd-Δμfg) (*liter atm*)

Calculate p from 
	p = 1- Δμ(298, 1)/Δv(298,1)

p=1-Δμ/Δv  (*atm*)

WorkBookT23.3 Calculate approximately the vapor pressure of water

The exact equation for the dependence of the boiling temperature on pressure is 
	(dp/dT) = ((Δh(T))/(TΔv(T)))						(1)
The system has only one variable and I take it to be T, the coexistence temperature. It is common (in textbooks) to assume that 
• Δh is independent of temperature.
• Δv = v(vapor) - v(liquid) ≃ v(vapor) ≃ R T/p 
• If these approximations are made Eq. 1 can be integrated to yield
	ln((p/(Subscript[p, i]))) = -  ((Δh)/(R )) ( (1/T) - (1/(Subscript[T, i])))				(2)
I want to test how well (or how poorly) this equation works.

Cell 1: Calculate the heat of evaporation Δh at 298 K

Off[Remove::rmnsm]
Remove["Global`*"] ;
Needs["Utilities`CleanSlate`"];
CleanSlate[];

The heat of formation of (Subscript[H, 2])O(g) under standard conditions is:

Δhfg=-241.8264 (*kJ/mol*);

The heat of formation of (Subscript[H, 2])O(l), under standard conditions, is:

Δhfl= -285.8304 (*kJ/mol*);

The heat of evaporation is the heat of formation of water vapor minus the heat of formation of water liquid (all at 298.15 K and p = 1 atm).

Δh=Δhfg-Δhfl (*kJ/mol*)

Cell 2: Calculate the vapor pressure

Equation (2) gives  
	p/(Subscript[p, i]) = Exp[- ((Δh)/R) ((1/T)-(1/(Subscript[T, i]))) ]   		(3)
The  vapor pressure of water at  100 degree Celsius  (which is 372.15 K) is 1 atm.  Therefore, I took 
	(Subscript[T, i]) =373.15 K and (Subscript[p, i]) = 1 atm.

Calculate the vapor pressure from Eq. 3. Since Δh is in kJ/mol use 
	R = 8.31451 (10^(-3)) kJ/mol K.

Clear[p]
p[T_]:=With[{R=8.31451 10^-3},
Exp[-(Δh/R)( 1/T -1/373.15)] ] (*dimensionless; p is divided by 1 atm*)

p[x]

Generate a list of temperatures

temp=Table[273.15+t,{t,10,100,10}]

Calculate, by using Eq. 3, the vapor pressures of water at these temperatures :

pCalc=p/@temp

Cell 3: Compare the calculated vapor pressure with the measured values

Compare these results with the experimental values (page D168 in CRC Handbook of Physics and Chemistry):

pExp={9.209,17.535,31.824,55.324,92.51,149.38,233.7,355.1,525.76,760.}/760 (*atm*);

error=(pExp-pCalc)*100/pCalc;

TableForm[Transpose[{temp,pExp,pCalc, error}],
TableHeadings->{None,{"T K", "(Subscript[p, exp])","(Subscript[p, calc])","%error"}}]

This is pretty bad agreement.  The approximations made do not work.

WorkBookT.23.4  Test the Clapeyron equation for (Subscript[NH, 3]) evaporation.

Off[Remove::rmnsm]
Remove["Global`*"]
Off[General::spell1]

• I will use data for ammonia, taken from an article by P.Davies on p. 33  in the book Thermodynamic Functions of Gases, vol. 1, edited by F. Din, published by Butterworth, London, 1962.  
• Davies collected a remarkably detailed set of data concerning ammonia liquid-vapor equilibrium: the molar volume of the liquid, the molar liquid of the gas, the heat of evaporation as a function of the coexistence temperature.  
• This allows we to verify how good the Clapeyron equation is. 
• Surprisingly, for such careful work, I found that one of the equations given by Davis is very erroneous. Luckily the experimental data was available and I made a correct fit.

Cell 1: The vapor pressure versus temperature

I use an expression for vapor pressure versus temperature given by Eq. 22 , p.47 in Din vol. 1.  This equation was obtained by fitting the  vapor pressure measured at different temperatures

Input the measured vapor pressure:

pVaporExp={0.08553,0.1756,0.3347,0.5982,1.0124,1.6320,2.528,3.769,5.446,7.652,10.485,14.063,18.49,23.90,30.42,38.19,47.40,58.145,70.60,84.94,101.5};

These vapor pressures correspond to the following temperatures

temp=Table[t,{t,200,400,10}] ;

This data was fitted by Davis with the following equation

p[T_]:=10^(27.376004-1914.9569/T -8.4598324 Log[10,T]+
2.393909*10^(-3) T+2.955214*10^(-6) T^2) (*vapor pressure in atm*)

I have the data and I have an equation representing it.  I will test how good the fit is.  Calculate the vapor pressure provided by the fit

pVaporCalc=p/@temp;

Make a table with the calculated pressure, the experimental one and the error of the fit.

toTable=MapThread[{#1,#2,#3,100(#2-#3)/#3}&,{temp,pVaporCalc,pVaporExp}];TableForm[toTable,TableHeadings->{None,{"T K", "pCalc","pExp","error"}}]

The equation for pressure (implemented by the function p[T]) is very good. Make a plot of p[T] versus T, which is the vapor-liquid coexistence line. Show the data (as dots) and the equation (as a line).

toPlot=MapThread[{#1,#2}&,{temp,pVaporExp}];
pl1=ListPlot[toPlot,
PlotStyle->{PointSize[0.02]},
DisplayFunction->Identity];
pl2=Plot[p[T],{T,200,400},
DisplayFunction->Identity];
Show[{pl1,pl2},DisplayFunction->$DisplayFunction,
AxesLabel->{"T K", "p atm"},
PlotLabel->"(Subscript[NH, 3]) vapor pressure"]

A given state {p, T} is a point in this graph.  If that point is at the left of the line (Subscript[NH, 3]) is liquid; if it is at the right it is solid.  If the point is on the line the system has two phases (liquid and vapor) in equilibrium with each other.

Cell 2: The molar volume of the gas coexisting with the liquid

Use equation 28 in Din, vol. 1, p. 52 which gives the dependence of the gas volume (in (cm^3)/mol) on temperature, when the two phases are in equilibrium.

vGasCoex[T_]:=Module[{a,b,c},
a=300(6.4644/T -0.106887+ 0.0356803 Log[10,T]) ;
b=0.0862366 Sqrt[406.1-T] ;
c= 0.002667*(406.1- T) ;
 17*10^(a+b+c) (*volume in (cm^3)/mol*)]

Display the equation:

vGasCoex[T]

Calculate the volume of the gas ((cm^3)/mol) on the coexistence curve (i.e. when the vapor is in equilibrium with the liquid) at the temperatures at which I have data:

vGasCoexCalc= (vGasCoex/@temp ) (*(cm^3)/mol*);

Compare to the experimental data (Din, vol. 1, p. 89):

vGasCoexExp={189730,96858,53070,30895,18905,12088,8029,5505,3882,2803,2066,1545,1174,900.1,696.1,540.4,419.3,323.9,247.0,182.6,124.5} (*(cm^3)/mol*);

Make a table showing the calculated volumes, the experimental one and the error of the equation.

Print["the molar volume is in (cm^3)/mole"]
TableForm[MapThread[{#1,#2,#3,100*(#3-#2)/#3}&,{temp,vGasCoexCalc,vGasCoexExp}],
TableHeadings->{None,{" T K","vCalc","vExp","error"}}]

The equation provided by Davis is fairly good.

I will make a temperature volume phase diagram. For this I need to plot the temperature versus volume of gas

pl1=ListPlot[ Transpose[{vGasCoexCalc,temp}], PlotStyle->PointSize[0.02],
DisplayFunction->Identity ];Show[pl1,AxesLabel->{"v (cm^3)/mol","T K"},PlotLabel->"volume of saturated vapor" ]

This graph may seem odd: the volume of the coexisting vapor decreases when the temperature increases. This is not possible for a gas in a container (no liquid).  In the case of coexistence the temperature is the independent variable.  At low T the vapor pressure is very low and therefore the molar volume is very high (as the graph shows)

Cell 3: The molar volume of the liquid coexisting with the gas

Din gives Eq. 27 at p. 51 if volume 1, which is supposed to give the molar volume of the liquid when it coexists with the vapor.

vLiqCoex[T_(*Kelvin*)]:=Module[{t,a,b,c},
t=T-273.15;
a=4.2830+0.81305*Sqrt[133-t]-0.0082861*(133-t);
b=1+0.424805*Sqrt[133-t]-0.015938*(133-t);
17*a/b (*(cm^3)/mol*)]

Display the equation:

vLiqCoex[T]

This is an odd fit and I will show that it does not work well.  It is not clear to me why a square root was used. The volume is a complex number if T  > 406.15.  It is not clear why 406.15 K is in the equation.

Calculate the volumes at the temperatures for which I have data:

vLiqCoexCalc=vLiqCoex/@temp;

The measured volumes, at the same temperatures are

vLiqCoexExp={23.37,23.74,24.13,24.55,24.99,25.46,25.96,26.50,27.07,27.70,28.39,29.14,29.98,30.92,31.99,33.24,34.73,36.58,39.03,42.61,49.45} (*(cm^3)/mol*);

Compare the calculated volumes with the measured ones:

TableForm[MapThread[{#1,#2,100(#1-#2)/#2}&,{vLiqCoexCalc,vLiqCoexExp}],
TableHeadings->{None,{"(Subscript[v, calc])","(Subscript[v, exp])","% error"}}]

The measured volumes are very different from the calculated ones! There is an error in the equation given by Davis/Din. I will use Mathematica function Fit to make my own fit of data. 
First arrange the data in the form required by the function Fit

toFit=Transpose[{temp,vLiqCoexExp}];

Fit the data with a  linear combination of the following terms:  a constant, 405.6 -T, (((405.6-T))^2), 1/(405.6 - T) and 1/(((405.6-T))^2)

theFit=Fit[toFit,{1,405.6-T,(405.6-T)^2,1/(405.5-T),1/(405.6-T)^2},T]

This is not trouble free because the function has a pole at T = 406.15 K.  However, the coexistence temperatures I examine here are smaller than 406 K, as you can see from the list  below.

temp

vLiqCalc[x_]:=theFit/.T->x

Check if it works:

vLiqCalc[z]

Use the function to calculate the volume of the liquid at coexistence, for the temperatures for which I have data. Compare to experiment by calculating the error.

fitList=vLiqCalc/@temp;
errorV=100*(vLiqCoexExp-fitList)/vLiqCoexExp;
TableForm[Transpose[{temp,vLiqCoexExp,vLiqCalc[temp],errorV}],
TableHeadings->{None,{"T K","(Subscript[V, exp]) (cm^3)/mol","(Subscript[V, calc]) (cm^3)/mol","%error"}}]

The fitting function is fairly good.  I am surprised that the fit is  good even at 400 K, which is very close to the fatal 405.6.

Make a plot of temperature versus the molar volume of the liquid coexisting with the vapor.

pl2=ListPlot[Transpose[{ vLiqCoexExp,temp}],AxesLabel->{"v (cm^3)/mol","T k"},
PlotStyle->PointSize[0.02]]

Cell 4: The heat of evaporation

The dependence of the heat of evaporation on temperature is given in Din vol. 1, p. 58, Eq. 33:

heat[T_]:=Module[{},
t=T-273.15;
0.239*17*(137.91*Sqrt[133-t]-2.466*(133-t))(*cal/mol*)]

heat[T]

This equation is valid for temperatures between - 42 C and 52 C  (231 K and 325K).

The experimental data for heat of evaporation is

hEvapExp={5982,5892,5795,5686,5567,5446,5318,5182,5037,4882,4716,4539,4345,4146,3916,3656,3366,3027,2609,2079,1269} (*cal/mol*);

Test the equation:

heatCalc=(heat/@temp);

error=100*(heatCalc-hEvapExp)/hEvapExp;

TableForm[Transpose[{temp,heatCalc, hEvapExp,error}],
TableHeadings->{None,{"T K","(Subscript[Δh, calc])","(Subscript[Δh, exp])","%error"}}]

The equation is accurate. The only bad point is T=400 K where the error is 4.6%

Cell 5: Test Clapeyron equation

I have data at the temperatures

temp=Table[T,{T,200,400,10}];

The equation below gives the vapor pressure in atm. I have already tested it.

p[T_]:=10^(27.376004-1914.9569/T -8.4598324 Log[10,T]+
2.393909*10^(-3) T+2.955214*10^(-6) T^2) (*vapor pressure in atm*)

The experimental molar volumes of the liquid at coexistence for the temperatures in temp are:

vLiqCoexExp={23.37,23.74,24.13,24.55,24.99,25.46,25.96,26.50,27.07,27.70,28.39,29.14,29.98,30.92,31.99,33.24,34.73,36.58,39.03,42.61,49.45} (*(cm^3)/mol*);

The molar volumes of the vapor at coexistence for the temperatures in the list temp are:

vGasCoexExp={189730,96858,53070,30895,18905,12088,8029,5505,3882,2803,2066,1545,1174,900.1,696.1,540.4,419.3,323.9,247.0,182.6,124.5} (*(cm^3)/mol*);

The measured heats of evaporation at the temperatures in the list temp are:

hEvapExp={5982,5892,5795,5686,5567,5446,5318,5182,5037,4882,4716,4539,4345,4146,3916,3656,3366,3027,2609,2079,1269} (*cal/mol*);

Test the equation 
	dp/dT =Δh/ TΔv			(1)

Every quantity entering in Eq. 1 has been measured.  Let us see if the measured values satisfy this equation:

• Calculate dp/dT from equation (1).  
• Δh(T) is given in cal/mol and I convert it to liter atm by using 
	1 cal = 0.041297 liter atm.  
• The volumes are given in (cm^3)/mol and I convert them in liter/mol (1 liter =(10^3)(cm^3)).  
• After these conversions  Δh/TΔv will have units of atm/K, which is what I want.

The right-hand side of the  Clapeyron equation 
	dp/dT =Δh/ TΔv 
is calculated below by using experimental data for Δh and Δv

dpdtClap=0.041297*hEvapExp/(temp*(vGasCoexExp-vLiqCoexExp)*10^-3) (*atm/K*);

I can also calculate dp/dT from the equation that fitted the measured vapor pressure to the coexistence temperature T:

dpdT=D[p[T],T]/.T->temp  (*atm/K*);

TableForm[MapThread[{#1,#2,#3,100*(#3-#2)/#3}&,{temp,dpdtClap,dpdT}],
TableHeadings->{None,{"T K","(dp/dT)calc","(dp/dT)exp","%error"}}]

Obviously the equation is very accurate.

WorkBookT.23.5 Integrate the Clausius-Clapeyron equation to determine the vapor pressure

I want to integrate the equation 
	(dp/dT) =((Δh(T))/(T Δv(T)))
and calculate the vapor pressure dependence on temperature.  This is a first order differential equation and to solve it I need to know the pressure at one temperature.   I know how Δh(T) and Δv(T)  depend on temperature.

Cell 1:  Collect that data  you need

The gas molar volume at coexistence is given by :

vGasCoex[T_]:=Module[{a,b,c},
a=300(6.4644/T -0.106887+ 0.0356803 Log[10,T]) ;
b=0.0862366 Sqrt[406.1-T] ;
c= 0.002667*(406.1- T) ;
 17*10^(a+b+c) (*volume in (cm^3)/mol*)]

vGasCoex[T]

The liquid molar volume at coexistence is:

vLiq[T_]:=35.38385  + 170.75158 /(405.5 - T) - 
  515.23394/(405.6 - T)^2 - 0.100695*
   (405.6 - T) + 0.00018915*(405.6 - T)^2

vLiq[x]

The change in molar volume caused by evaporation is:

Δv[T_]:=vGasCoex[T]-vLiq[T] (*(cm^3)/mol*)

The heat of vaporization is given by:

Clear[Δh]
Δh[T_]:= Module[{t},
t=T-273.15;
0.239*17*(137.91*Sqrt[133-t]-2.466*(133-t))(*cal/mol*)]

Cell 2:  Solve the differential equation

I have 
	(dp/dT) =((Δh(T))/(T Δv(T)))
Regard this as a differential equation  which is to be solved to obtain p[T]

Write down the differential equation to be solved numerically:

eq={pre'[T]==0.041297*Δh[T]/(T*Δv[T]*10^-3),pre[200]==0.08553};

• pre[T] is the function that gives the vapor pressure versus T. 
• Because the heat of evaporation is given by the function Δh(T) in calories I use 0.041297 to convert it to liter atm. 
• Δv is given in (cm^3)/mol and I convert it to liter ((10^(-3))liter=(cm^3)). 
• I use the initial condition:  p = 0.08553 atm for T = 200 K, which I know from experiment

Solve the equation numerically

sol=NDSolve[eq,pre,{T,200,400}]//Flatten

Turn the interpolating function into an usable function pC  (pressure by Clapeyron)

pC=pre/.sol

Calculate the vapor pressure given by solving the differential equation at the temperatures for which I have experimental data:

pCList=pC/@temp

Cell 3:  Compare to experiment

Compare to the experimental results at the same temperatures:

pVaporExp={0.08553,0.1756,0.3347,0.5982,1.0124,1.6320,2.528,3.769,5.446,7.652,10.485,14.063,18.49,23.90,30.42,38.19,47.40,58.145,70.60,84.94,101.5};

TableForm[MapThread[{#1,#2,#3,100*(#2-#3)/#3}&,{temp,pCList,pVaporExp}] ,
TableHeadings->{None,{"T K","(Subscript[p, calc]) atm", "\!\(\*SubscriptBox[\(p\), \(\(exp\)\(\\ \)\)]\)atm", "%error"}}]

This is an excellent agreement, except at the highest temperature where I know that Δh is not accurately rendered by the equation Δh(T).