WorkBookT3. Equations of State  © 2022 Horia Metiu

Remove["Global`*"];
Needs["Utilities`CleanSlate`"]
CleanSlate[]

Cell 1: Calculate the pressure when you know v and T.

Cell 1.1:  Explain the problem and the plan for solution

• Use the ideal gas, the van der Waals, the Beattie-Bridgeman, and the Benedict-Webb-Rubin equations to calculate the pressure of ethane at the temperature 
	T = 350 K, 
when the molar volume has the values:
 	v = {28.588, 5.6119, 2.7387, 0.1372, 0.0984, 0.0868,  0.30716, 0.07089, 0.07022} liter/mole. 
• Compare the results of these calculations with the experimentally measured pressures:
	 p =  {1, 5, 10, 100, 150, 200, 450, 475, 500} atm  (see Appendix 2)
Note: You read the data in these two lists as follows:  when the volume is 5.6119 liter/mole the pressure is p = 5 atm, etc.

I know T and v and I want to calculate p. All I have to do is to use the appropriate values for T and v in the equation of state p = f(T, v).  For example, if   					
	v = 28.588 liter/mol and T = 350 K, 
then the pressure is 
	p = f(350, 28.588). 
I am required to do this calculation for four different equations of state. The ideal gas law (Eq. 2.2 in textbook)
	p = f(T, v) = R T/v; 									(1)
the van der Waals equation  (Eq. 2.9)
	p = f(T, v) = (RT/(v-b))-(a/(v^2)); 									(2)
the Beattie-Bridgeman equation (Eq. 2.10)
	 p = f(T, v) = ((R T)/(v^2)) (1 -  (c/(v (T^3)))) (v +B-((b B)/v)) -   (A/(v^2))(1-(a/v))  			(3)
and the Benedict-Webb-Rubin equation  (Eq. 2.11): 
	p = R T ρ + (B R T-A-(ρ^2)(Subscript[C, 0])T-2  + (ρ^3)(b R T-a)+ a α (ρ^6)  + 
	(ρ^3)cT-2 (1 + γ (ρ^2)) exp[-γ (ρ^2)]							(4)
 where  ρ =1/v.  The constants appearing in these equations are given in Appendices  4-6 of the textbook.

 I break up the calculation into pieces. 
• Define functions that calculate the pressure according to the four equations given above.  
• Get the constants I need from the Appendices.  
• Input the volume and pressure data.
• Calculate the pressures corresponding to the volumes given in the problem (the temperature is fixed to 350 K).   
• Compare the pressures calculated with the equations of state to the measured pressures.
• Examine why the van der Waals equation fails so miserably at high pressure.

Cell 1.2: Define the functions that calculate the pressure

I define the functions 
	pIdeal[T, v], 
	pvdw[T, v, parameters],
	pBB[T, v , parameters] , 
	pBWR[T, v, parameters]
that give the pressure when I know v and T.  
	pIdeal uses the ideal gas law, 
	pvdw uses  the van der Waals equation, 
	pBB uses the Beattie-Bridgeman equation, 
	pBWR uses the Benedict-Webb-Rubin equation.

Define the function pIdeal[T,v] that calculates the pressure from the ideal gas law 
	p = R T/v

Clear[pIdeal]
pIdeal[T_(*Kelvin*),v_(*liters*)]:=With[{R=0.08205(*liter atm/mol*)},
		p= R T/v]

Test that it works

{pIdeal[T,v],pIdeal[278,24]}

Mathematica note:  I use function With provided by Mathematica because it will keep the quantity R = 0.08205 private ( known only in the expression inside With).

Define the function 
	pvdw[T,v, parameters]  
that calculates the pressure from the van der Waals equation 
	p = f(T, v) = (RT/(v-b))-(a/(v^2))
I define this function so that the parameters are given through replacement rules such as  
	parameters = {a→ number1, b→number2}

Clear[pvdw]
pvdw[T_,v_,parameters_]:=Module[{R=0.08205 (*liter atm/mol*) },
   Clear[a,b];
		(R T/(v-b)-a/v^2)/.parameters];
SetAttributes[pvdw,HoldAll];

The attribute HoldAll prevents the argument of the function from being evaluated.  Test that the equation works well.  The parameters for ethane are (see Appendix 4)
	a = 5.489 (liter^2) atm/mole
	b = 0.0638 liter/mole

parametersvdw={a->5.489 ,b->0.0638}

{pvdw[T,v,{}],pvdw[T,v,parametersvdw]}

Below I calculate the pressure of ethane at T= 298 K and v = 24 liter . If the function is correct this should be close to 1 atm.

pvdw[298,24,{a->5.489 ,b->0.0638}]

Define the function pBB[T, v ,parameters] to calculate the pressure from the Beattie-Bridgeman equation, Eq. 3: 
 	p = f(T, v) = ((R T)/(v^2)) (1 -  (c/(v (T^3)))) (v +B-((b B)/v)) -   (A/(v^2))(1-(a/v))

Clear[pBB]
pBB[T_(*liter/mol*),v_(*liter/mol*),parameters_]:=Module[
{ R=0.08205(*liter atm/mol K*),p1,p2,p3},
     Clear[A,a,B,b,c] ;  
    p1=1-c/(v T^3);
    p2=v+B-b B/v;
    p3=A*(1-a/v);
   ( (R T p1 p2-p3)/v^2)/.parameters (*p in atm*)];
SetAttributes[pBB,HoldAll];

The constants in the Beattie-Bridgeman equation are given through replacement rules

parametersBB={A->5.8800,B->0.094,a->0.05861,b->0.01915,c->90.0 10^4}

At room temperature and molar volume of 24 liter the pressure should be close to 1 atm. Verify:

pBB[298,24,parametersBB]

Define a function 
	pBWR[T, v, parameters] 
which calculates the pressure by using the Benedict-Web-Rubin equation (Eq. 4):
	p = R T ρ + ( B R T-A-(Subscript[C, 0])T-2) ρ2 +(b R T-a) ρ3+a α ρ6 + 
	+ c (ρ^3) (T^(-2)) (1+γ ρ2) exp[-γ (ρ^2)] 
Here ρ = 1/v where v is the molar volume.

Clear[pBWR]
pBWR[T_,v_,parameters_]:=Module[
{R=0.08205,p1,p2,p3,p4,p5,ρ},
Clear[A,B,C0,a,b,c,α,γ];
 ρ=1/v;
p1=R T ρ;
p2=(B R T-A-C0 T^-2)ρ^2;
p3=(b R T-a)ρ^3;
p4=a α ρ^6;
p5= c T^-2 ρ^3 (1+γ ρ^2) Exp[-γ ρ^2];
(p1+p2+p3+p4+p5)/.parameters];
SetAttributes[pBWR,HoldAll];

The parameters in the Benedict-Webb-Rubin equation, for ethane, are given in Appendix 6.

Clear[A,B,C0,a,b,c,γ]
parametersBWR={A->4.15556 ,B->6.27724 10^-2 ,C0->1.79592 10^5 , a-> 0.34516,b->1.1122 10^-2 ,c->3.2767 10^4 ,α-> 2.43389 10^-4,γ->1.18 10^-2 };

pBWR[298,24,parametersBWR]

Cell 1.3: Make lists with the experimental volume and pressure.

The molar volumes given in the list vExp, defined below correspond to the pressures given in the list  pExp. For example, if p = 10 atm, then v = 2.7387 liter/mol.  All data is for T = 350 K.

Clear[vExp]
vExp=  {28.588,5.6119,2.7387,0.1372,0.0984,0.0868,0.0716,0.07089,0.07022} (* liter/mol*);
Clear[pExp]
pExp={1, 5, 10,100,150,200,450,475,500} (* atm*);

Cell 1.4: Calculate the pressures by using all four equations of state

To do this, I Map the appropriate equation of state into the list containing the values of the molar volumes.  The expression  
	pIdeal[350, #]&/@vExp 
calculates the pressure given by the ideal gas equation for the volumes in the list vExp, at a temperature of 350 K.  The expression  
	pvdw[350,#]&/@vExp 
does the same thing for the van der Waals equation, etc

Mathematica explanation:  f/@{a,b,c} maps the function f onto the list {a,b,c}:

f/@{a,b,c}

The syntax pIdeal[350, #]&/@{a,b,c} works as follow:

pIdeal[350,#]&/@{a,b,c}

Here  28.7 = R 350.  The function pIdeal[350,#]& is Mapped by the symbol /@ into the list  {a, b, c}.

Below, the list pIdealList contains the pressures calculated with the ideal gas law, pvdwList those calculated with the van der Waals equation, etc.

pIdealList=pIdeal[350,#]&/@vExp; 
pvdwList=pvdw[350,#,parametersvdw]&/@vExp;
pBBList=pBB[350,#,parametersBB]&/@vExp;
pBWRList=pBWR[350,#,parametersBWR]&/@vExp;

Cell 1.5: Make a table with the results of the calculations

I want to display in a table, the pressures I have calculated.  The information is in the lists 
	vExp, pExp, pIdealList, pvdwList, pBBList, pBWR.
I need to organize them so that I can display them with TableForm. The Transpose function does that:

toTable=Transpose[{vExp,pExp,pIdealList,pvdwList,pBBList,pBWRList}];Print["The temperature is T = 350K"];TableForm[toTable,TableHeadings->{None,{"(Subscript[v, exp])","(Subscript[p, exp])","(Subscript[p, ideal])",
			"(Subscript[p, vdW])","(Subscript[p, BB])","(Subscript[p, BWR])"}}]

Obviously by the time p = 100 atm the ideal gas law is completely wrong. The vdw equation is still OK.  At 500 atm only the Benedict-Webb-Rubin is correct. van der Waals is very bad if   p >100 atm.

Mathematica explanation: see below how the Transpose and TableForm work together. If the lists are long the syntax shown below won’t work

TableForm[{{a,b,c,d,e,f},{1,2,3,4,5,6}}]

Transpose[{{a,b,c,d,e,f},{1,2,3,4,5,6}}]

TableForm[Transpose[{{a,b,c,d,e,f},{1,2,3,4,5,6}}]]

The last syntax  is better for what I want.

Cell 1.6:  Calculate the error made by these calculations

Mathematica explanation: I use here the properties of operations with lists exemplified below

({a,b,c}-{A,B,C})/{x,y,z}

I calculate the errors in the pressures calculated by each equation of state

errorIdeal= (pExp-pIdealList)*100/pExp;  
errorvdw=(pExp-pvdwList)*100/pExp;
errorBB=(pExp-pBBList)*100/pExp;
errorBWR=(pExp-pBWRList)*100/pExp;

Display these results as a table:

toTable=Transpose[{pExp,errorIdeal,errorvdw,errorBB,errorBWR}];Print["The temperature is T = 350 K"];
		PaddedForm[TableForm[toTable,TableHeadings->{None,{"pExp","%errorId","%errorvdw","%errorBB","%errorBWR"}}],{6,3}]

PaddedForm controls the output (I don’t want more than 3 decimal accuracy). The errors are small at low pressure. The ideal gas law starts having trouble at 5 atm. The van der Waals equation is OK up to 100 atm but then breaks down severely.  Only the Benedict-Webb-Rubin equation is adequate through the whole pressure range.

Cell 1.7: Why is the van der Waals equation so bad when v is small?

When the molar volume v is close to the parameter b, the van der Waals equation gives rather poor results. To see why, look at the graph below.

Plot[pvdw[350,v,parametersvdw],{v,0.05,0.09},
	PlotLabel-> "         van der Waals equation" ,
AxesLabel->{"v liter","p atm"}]

When v < b, the pressure is negative!  It goes to infinity when v = b. This is absurd and it happens because the equation is not valid for these small molar volumes.  p =(RT/(v - b))- (a/(v^2))   cannot give good results even when v > b and v is  close to b  and it is completely wrong when  v < b.

Cell 2: Calculate the molar volume when you know p and T.

Cell 2.1: Formulate the problem and the plan for solving it.

Problem:  Calculate the molar volume of ethane, at a temperature of 400 K and the pressures 
 	{1, 2, 3, 5, 6, 7, 8, 9, 10} atm. 
 Use the van der Waals equation, and the Beattie-Bridgeman equation. The parameters for the van der Waals and the Beattie-Bridgeman equations are given in Appendices 4 and 5.

• When I know the pressure p and the temperature T, I can calculate the molar volume v from the equation of state  by solving p = f(T, v) for v.				
• I will use the van der Waals equation 
	p = R T/(v-b) - a/(v^2) 									(2)
• I will also use the Beattie-Bridgeman equation 
	 p = f(T, v) = ((R T)/(v^2)) (1 -  (c/(v (T^3)))) (v + B-((b B)/v)) -   (A/(v^2))(1 - (a/v))			(3)
• The constants appearing in these equations are given in Appendices 4 and 5. 
• The units are such that if I use the temperature in Kelvin, the pressure in atm and  R = 0.082 liter atm/mole K I get the volume in liter.
• Both these equations of state are polynomials is v.  Solving for v gives multiple solutions. One must decide, based on physical arguments,  which one is the molar volume of the gas  (see Ch. 3, §6-7). 

Mathematica strategy: 
• Define a procedure that calculates v from the equation p = f(T, v).  
• Calculate the molar volumes for the pressures indicated in the problem, by using all three equations of state. 
• Pick the solution that is physically reasonable.
• Display the results.

Cell 2.2: Define the van der Waals and Beattie-Bridgman equations

Clear[pvdw]
pvdw[T_,v_,parameters_]:=Module[{R=0.08205 (*liter atm/mol*) },
   Clear[a,b];
		(R T/(v-b)-a/v^2)/.parameters];
SetAttributes[pvdw,HoldAll];

The attribute HoldAll prevents the argument of the function from being evaluated.  Test that the equation works well.  The parameters for ethane are (see Appendix 4)
	a = 5.489 (liter^2) atm/mole
	b = 0.0638 liter/mole

parametersvdw={a->5.489 ,b->0.0638}

{pvdw[T,v,{}],pvdw[T,v,parametersvdw]}

Below I calculate the pressure of ethane at T= 298 K and v = 24 liter . If the function is correct this should be close to 1 atm.

pvdw[298,24,{a->5.489 ,b->0.0638}]

Define the function pBB[T, v , parameters] to calculate the pressure from the Beattie-Bridgeman equation, Eq. 3 
 	p = f(T, v) = ((R T)/(v^2)) (1 -  (c/(v (T^3)))) (v +B-((b B)/v)) -   (A/(v^2))(1-(a/v))

Clear[pBB]
pBB[T_(*K*),v_(*liter/mol*),parameters_]:=Module[
{ R=0.08205(*liter atm/mol K*),p1,p2,p3},
     Clear[A,a,B,b,c] ;  
    p1=1-c/(v T^3);
    p2=v+B-b B/v;
    p3=A*(1-a/v);
   ( (R T p1 p2-p3)/v^2)/.parameters (*p in atm*)];
SetAttributes[pBB,HoldAll];

The constants in the Beattie-Bridgeman equation, for ethane, are given below through replacement rules

parametersBB={A->5.8800,B->0.094,a->0.05861,b->0.01915,c->90.0 10^4}

{pBB[T,v,{}],pBB[T,v,parametersBB],pBB[298,24,parametersBB]}//TableForm

Cell 2.3: Define a function that solves p = f(T, v) for v.

If I know p and T, I can calculate v by solving the equation p = f(T, v). For example, for the van der Waals equation I must solve
		p = R T/(v-b) - a/(v^2)
for v. The function is a  cubic polynomial in v; solving for v will give three solutions.   I design a function that calculates v for the gas.

volumeG[T_,p_,f_,parameters_]:=
	Module[{R=0.082 (*liter atm/mol*),eq,sol},
Clear[a,b,c,A,B,C0,γ];
	eq=p==f[T,v ,parameters];
		v/.FindRoot[eq,{v,R T/p}] ];

• Mathematica explanation:  
• Mathematica has several functions that solve algebraic equations. 
• I prefer to use here FindRoot (To learn more about FindRoot see Help or Wolfram book, or the Mathematica lecture Algebra.nb posted on the web ).  The syntax for calling FindRoot is: 
	FindRoot[p == f[T, v], {v, (Subscript[v, 0])}] 
• To use the function FindRoot one must give p and T  numerical values.  
• The curly bracket {v, (Subscript[v, 0])} indicates that I want to solve for v and that a reasonable guess for the solution is  (Subscript[v, 0]) .  
• I use the ideal gas law 
	(Subscript[v, 0]) = R T/p
to get (Subscript[v, 0])
• The procedure used by FindRoot returns the solution "closest" to my guess, which is the molar volume of the gas.  In this way I avoid the problem of finding multiple solutions.  This procedure may have troubles near the critical point. 
• Exercise:  Figure out how you can use Solve to calculate v. Hint:  Use Solve after all quantities (R, T, p and the parameters) were given numerical values.  Otherwise the output is very long and incomprehensible.

Cell 2.4: The molar volume of ethane at 400 K  and  p = {1, 2, . . , 10} atm.

Generate a list of the pressure you want to calculate the molar volume for

p=Table[i,{i,1,10}]

Below, I map the function 
	volumeG[400, #, eqOfState, parameters]& 
into the list of pressures p.

The temperature is 400 K. In the first line below I use the van der Waals equation of state pvdw

volvdw=volumeG[400,#,pvdw,parametersvdw]&/@p

Use the Beattie - Bridgeman equation with appropriate parameters:

volBB=volumeG[400,#,pBB,parametersBB]&/@p

Cell 2.5: Display the results in a table.

I want to compare to the experimental results which are given in Appendix 2 and are:

vExp={32.726,16.315,10.845,8.1101,6.4687,5.3747,4.5933,4.0073,3.5513,3.1866};

Make a Table to display the results. Use the procedure  shown in WorkBook 3.1.  I will use this table making syntax throughout the lectures.

Print[ " T = 400 K  v is in  liter/mol" ] 
TableForm[Transpose[{p,volvdw,volBB, vExp}],
	TableHeadings->{None,{"p atm","(Subscript[v, vdw]) ","(Subscript[v, BB]) ","(Subscript[v, Exp])"}}]

The results for all equations, at these pressures and temperature, are very good

Cell 2.6:  Multiple solutions for the molar volume

In all realistic equations of state the pressure dependence on the molar volume v is a polynomial. Therefore, I expect the equation p = f(T, v) to give more than one root when solved for v at fixed T and p. This raises two questions:
 • which root gives the correct molar volume?
 • what is the meaning of multiple roots?

In this cell I examine these questions for the van der Waals equation. It turns out that
 • At high temperature and low pressure the equation has only one real root and this is the molar volume of the gas.
 • At lower temperatures and or high pressure the equation has three real roots. One corresponds to the liquid, one to the gas, and one for an unstable system.

Define again the van der Waals equation

Clear[pvdw]
pvdw[T_,v_,parameters_]:=Module[{R=0.08205 (*liter atm/mol*) },
   Clear[a,b];
		(R T/(v-b)-a/v^2)/.parameters];
SetAttributes[pvdw,HoldAll];

The parameters for ethane are (see Appendix 4)
	a = 5.489 (liter^2) atm/mole
	b = 0.0638 liter/mole

parametersvdw={a->5.489 ,b->0.0638}

Examine first the high temperature case.  Take p = 10 atm and T = 400 K. 
 I  solve p = (RT/(v-b))-(a/((v^2) )) for v )

solvdw=NSolve[400==pvdw[10,v,parametersvdw],v ]

• The equation is a third order polynomial in v and must have three solutions. 
• They must be either 
	• one real and two complex (they are complex conjugate)
	• three real 
• For p = 400 atm and T = 10 K  I get two complex roots and a real one.  
• The real solution is the molar volume of ethane when 
	p = 10 atm and T = 400 k.

Let’s plot p versus v at 400 K.  This is called an isotherm because T is constant.

Plot[pvdw[400,v,parametersvdw],{v,5,30},PlotLabel->"       the pressure dependence of p on v for T= 400K",AxesLabel->{"v liters","p atm"}]

Nothing surprising here: at high pressure the molar volume is small. You lower the pressure the molar volume in creases. To any given pressure corresponds on molar volume.

Now let’s use the van der Waals equation to plot the pressure as a function of the molar volume for a lower temperature isotherm.  Try T = 270 K.

Plot[pvdw[270,v,parametersvdw],{v,0.1,1}, PlotRange->All,AxesLabel->{"v liter","p atm"},PlotLabel-> "van der Waals isotherm at T=270K" ]

What is going on?  At the low temperature  I used here the isotherms look very different from what I had at high temperature. In particular, note that according to the graph for  pressures below ~30 atm three volumes are possible.

solvdw=NSolve[20==pvdw[270,v,parametersvdw],v ]

The van der Waals equation is telling us that something qualitatively different happens when the temperature is sufficiently low and the pressure is sufficiently high. The only explanation that any one could come up with is that the high molar volume corresponds to a gaseous ethane and the smallest molar volume corresponds to liquid ethane.  But, how can one interpret the middle values?

To make sense of the middle value let’s calculate the derivative of the pressure with the volume at 270 K and p = 25 atm

D[pvdw[270,v,parametersvdw],v]/.NSolve[25==pvdw[270,v,parametersvdw],v ,Reals]

I find that for the intermediate volume of  v = 0.175 (the middle value of v) 
	 ((∂p)/(∂v)) = 185.288 atm/liter
But this is absurd! If such a system existed it would have the strange property that it molar volume increases when the  pressure increases.  This has never been observed and because of this the middle value of the volume has been ruled unphysical.

Calculate the volume from the Beattie-Bridgeman equation.  I ask for positive, real roots.

NSolve[25==pBB[270,v,parametersBB]&&v>0,v] //Flatten

This is an equation of fourth order and has three positive solutions .  Look at the derivatives

D[pBB[270,v,parametersBB],v]/.NSolve[25==pBB[270,v,parametersBB]&&v>0,v]

• We see the same behavior as in the case of the van der Waals equation. 
• There are three roots.  The small one is the molar volume of the liquid, the largest one is the molar volume of the gas and the middle one corresponds to an unphysical situation (it can never be observed).   
• The plot of pressure versus volume (at T = 270 K) is similar to that obtained for the van der Waals equation.

Plot[pBB[270,v,parametersBB],{v,0.01,3},PlotRange->{-35,40},PlotLabel->"         Isotherm at T=270 K for Beattie-Bridgeman equation",AxesLabel->{"v liter/mole","p atm"}]

Cell 3:  Calculate the temperature, when you know the pressure and molar volume

Problem:  You have a gas at equilibrium whose molar volume is 23.595  liter/atm. At what temperatures is the pressure equal to  p = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10} atm?  Use the Beattie-Bridgeman equation of state (Eq. 2.10).

Plan:  I know v and p, and I calculate T by solving the equation p = f(T, v), where f is the equation of state.  
Mathematica strategy: This parallels the calculation of v in WorkBook 3.2. But here I calculate T instead of v.

Cell 3.1: The parameters in the Beattie-Bridgeman equation for ethane

The parameters in the Beattie-Bridgeman equation for ethane are given in Appendix  5.  They are

parametersBB={A->5.88,a->0.05861,B->0.094,
 b->0.01915, c->90.00*10^4};

Cell 3.2: Calculate the temperature  when you know v and p

Now I need a function 
	t[p, v, f, parameters] 
that calculates the temperature by solving the equation 
	p ==  f(T,v)
where f is the equation of state and p and v have known numerical values. In what follows I fix the molar volume to v = 23.595 liter/mol. This means that the number of moles/volume does not change (no gas goes in or out of the chamber and the volume of the chamber does not change)

If you have run the previous cells the Beattie-Bridgeman equation of state is in the memory.  Verify that

?pBB

Clear[t]
t[v_,p_,f_, parameters_]:=Module[{R=0.082 (*liter atm/mol degK*),sol},Clear[a,b,c,A,B];
T/.FindRoot[p==f[T,v,parameters],{T,p v/R}]
			];

t[23.595,1,pBB,parametersBB]

Cell 3.3:  At what temperatures is v = 23.595 mole/liter and the pressures p= {1, 2, ....,10} atmosphere

Perform a calculation of the temperature corresponding to the pressures 
	p = {1, 2, . . ., 10}

p=Table[i,{i,1,10}];

Calculate the temperature by mapping the function t in the list p

temp=t[23.595(*v in liter/mole*),#,pBB,parametersBB]&/@p;

Make a table with the results

TableForm[Transpose[{temp,p}],
	TableHeadings->{None,{"T K","p atm"}}]

Note that I have to increase the temperature quite a bit to get the pressure up by one atm.

ListPlot[Transpose[{temp,p}],
PlotStyle->PointSize[0.02],
AxesLabel->{"temperature K","pressure atm"}]

Should you be surprised that the dependence is a straight line?