© 2022 Horia Metiu

WorkBookT.S3. Make your own equation of state

WorkBookT.S3.1 A crude way of getting the virial coefficients

Cell 1: Define the problem and the data.

This solves the problem posed in Complement 3.1 by the method explained at §13.  
• At a temperature of 350 K and the pressures
	 p = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10} atm, 				(1)
the molar volumes of ethane are (see Appendix 2). The volume in liter/atm is
	 v = {28.588, 14.229, 9.441,7.047 ,5.612, 4.655, 3.970, 3.457, 3.058, 2.739}   (2)
• Determine the virial coefficients (Subscript[B, 2]), (Subscript[B, 3]), and (Subscript[B, 4]) in the virial equation 
 	p  = R T{ρ + (Subscript[B, 2]) (ρ^2) + (Subscript[B, 3]) (ρ^3) + (Subscript[B, 4]) (ρ^4)}				(3)
by using the "simple and unsatisfactory procedure" described in Complement 3.1, §13.  Here ρ = 1/v.

Plan for solving the problem: 
•  I have three unknown parameters, (Subscript[B, 2]), (Subscript[B, 3]) and (Subscript[B, 4]).  I pick up three pairs of data points
 	 	{p = 2, atm ρ = 1/14.229 mol/liter}, 		
  		{p = 4 atm, ρ = 1/7.0747 mol/liter} 
  and 
  		 {p = 6 atm, ρ = 1/4.665   mol/liter}.  
  • For each pair {(Subscript[p, i,]) (Subscript[ρ, i])}  I write the equation
		(Subscript[p, i])  = RT{(Subscript[ρ, i]) + (Subscript[B, 2]) ((Subscript[ρ, i])^2) + (Subscript[B, 3])  ((Subscript[ρ, i])^3) + (Subscript[B, 4]) ((Subscript[ρ, i])^4)}			(4)
Having three {p, ρ} pairs gives me three equations for (Subscript[B, 2]), (Subscript[B, 3]) and (Subscript[B, 4])
• I use function Solve to solve these equations for (Subscript[B, 2]), (Subscript[B, 3]) and (Subscript[B, 4]) 
• I use equation (4), with the coefficients (Subscript[B, 2]), (Subscript[B, 3]) and (Subscript[B, 4]) determined above, to calculate the pressures corresponding to the molar volumes given in (2) and compare the calculated pressures to the experimental values given in (1). 
• Calculate the errors made by the equation (determine how good the fit is)

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

Cell 2.  Make lists with the data

Create two lists, one with  the pressures for which I have data

dataP=Range[1,10](*atm*);

and one for the corresponding molar volumes (measured at T = 350K) taken from Appendix 2 (for ethane)

datav={28.588,14.229,9.441,7.047,5.612,4.655,3.970,3.457,3.058,2.739}(*liter/mole*);

The molar volume is 28.588 liter/mol when the pressure is 1 atm, etc.

I need only three data points to determine three coefficients.  I decided to take the second, the fourth, and the sixth data point:

shortDatav={datav[[2]],datav[[4]],datav[[6]]};

shortDataP={dataP[[2]],dataP[[4]],dataP[[6]] };

I want to work with the molar density ρ.  Calculate the short list of densities

shortDataρ=1/shortDatav

Cell 3:  Write the equations for B2, B3, and  B4, and solve them.

I want to generate the equations
	(Subscript[p, i]) = RT {(Subscript[ρ, i]) + (Subscript[B, 2]) (Subscript[ρ, i]^2)+ (Subscript[B, 3]) (Subscript[ρ, i]^3)+ (Subscript[B, 4]) (Subscript[ρ, i]^4)}		
where i = 1, 2, or 3, the quantities (Subscript[p, 1]), (Subscript[p, 2])and (Subscript[p, 3]) are three pressures in the list
 	shortDataP = {2, 4, 6}, 
 and  (Subscript[ρ, 1]),(Subscript[ρ, 2]) and (Subscript[ρ, 3]) are the three densities in the list 
	shortDataρ ={0.070279, 0.141904,  0.214823 }.

eq=.
eq= With[{R=0.082067  (*liter atm/mol K*),T=350},(#1==R T(#2 + B[2] #2^2 +B[3] #2^3+ B[4] #2^4))]&

Check that this function generates the equation I want when applied to a pair {p, ρ}

MapThread[eq,{{p},{ρ}}]

Mathematica explanation: the expression below uses MapThread to replace #1 in eq with a pressure and #2 with the corresponding molar density

eqs=.
eqs=MapThread[eq,{shortDataP,shortDataρ}]

Solve the system of equations for B1, B2, and B3.

sol=Solve[eqs,{B[2],B[3],B[4]}]//Flatten

Cell 4:  Define the virial equation for pressure

Define the function that gives the pressure according to virial equation
		p = RT {ρ + (Subscript[B, 2]) (ρ^2) + (Subscript[B, 3]) (ρ^3) + (Subscript[B, 4]) (ρ^4)}

p=.
p[ρ_]:=With[{R=0.082067 (*liter atm/K mole*),T=350 (*K*)},
R T (ρ +B[2] ρ^2 +B[3] ρ^3 +B[4] ρ^4)/.sol
(*gives the pressure as a function of ρ for T=350 K *)]

Check that it works properly

p[ρ]

Cell 5:  Compare the calculated pressures to the experimental ones

Make a list of molar densities ρ = 1/v:

dataρ=1/datav

Calculate the pressure corresponding to the densities in the list dataρ:

calcP=p/@dataρ

The correct result is reproduced below

dataP

Let us calculate the error made by calculating the pressures with the virial equation in which the coefficients were obtained with the naive method using three data points. The percent error is defined as
 	((Subscript[p, calc])-(Subscript[p, exp]))×100/(Subscript[p, exp])
Below we use MapThread  with an anonymous function to create a list in which each element is {calculated pressure, measured pressure, percent error}.

list=MapThread[{#1,#2,(#2-#1)*100/#2}&,{calcP,dataP}];

TableForm[list,TableHeadings->{None,{"P calc","P data","%error"}}]

We got perfect results for the three states used in fitting, and we are not doing badly for the other points.

WorkBookT.C3.2 Least squares fitting to get the virial coefficients

I have the following data for the pressure and the molar volume of ethane at T = 350 K (from Appendix 2):
{p, v} = {{1, 28.588}, {2, 14.229},{ 3, 9.441}, {4, 7.0477}, {5, 5.6119}, {6 , 4.6547}, {7, 3.9704}, {8, 3.457}, {9, 3.0581}, {10, 2.738}           				(1)
The units are atm and liter/mol, respectively.  I read the list as follows:  if p = 2 atm then v =14.229 liter/mol, etc. Use the least-square fitting method and this data to find the coefficients (Subscript[B, 2]), (Subscript[B, 3]), and (Subscript[B, 4]) in the virial equation 
	p = R T {ρ + (Subscript[B, 2 ])(ρ^2) + (Subscript[B, 3])(ρ^3)+(Subscript[B, 4])(ρ^(4 ))}					(2)

Plan for solving the problem: The least squares fitting method is described in Complement 3.1 in §14-§16.  I review it here.  If (Subscript[p, i]) is the experimental pressure when the molar volume is (Subscript[v, i]) and the temperature is T = 350K, then the global error is
	e((Subscript[B, 1]), (Subscript[B, 2]), (Subscript[B, 3])) =(∑_(i=1)^n)(([(Subscript[p, i])-f(T=350,(Subscript[v, i]))])^2)      
where f(T,v) is the virial equation of state
	 p = f(T,v) = R T(ρ  +(Subscript[B, 2])(ρ^2)+(Subscript[B, 3])(ρ^3)+(Subscript[B, 4])(ρ^4)).  
Here n is the number of data points. The best values for the virial coefficients are those that minimize the error e((Subscript[B, 1]), (Subscript[B, 2]), (Subscript[B, 3])) .  So I need to find the minimum of this function.

Mathematica strategy:  
•  Create lists with the data.  
•  Create a function that gives the pressure according to the virial equation.  
• Define a function that gives the global error.  
•  Minimize this function with respect to B[2], B[3], B[4], by using FindMinimum. 
• Test the virial equation with the coefficients determined above, by using it to calculate the pressures (for the given volumes) and comparing the results to the measurements. • • Make a table to display the results.

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

Cell 1: Create lists with the data.

The list of pressures:

pList=Range[10];

The list of molar volumes:

vList={28.588,14.228,9.441,7.0477,5.6119,4.6547,3.9704,3.457,3.0581,2.738};

The list of molar densities

ρList=1/vList;

Cell 2:  Create a function corresponding to the virial equation.

We need to create the function 
	pVirialρ[ρ, T, B2, B3, B4],   
that returns 
         R T { ρ + (Subscript[B2,  ])(ρ^2) + B3  (ρ^3)+B4 (ρ^(4 ))}

Clear[pVirial]
pVirial[T_,ρ_, B2_,B3_,B4_]:=Module[{R=0.082067 (*liter atm/mol K*)},
	Clear[B1,B2,B3];
	R T (ρ +B2 ρ^2 +B3 ρ^3 +B4 ρ^4) ];
	SetAttributes[pVirial,HoldAll]

Check out that pVirial forms the correct equation:

Clear[T,ρ, A2, A3, A4]
pVirial[T,ρ,A2,A3,A4]

The function pVirial  does return what I want.

Cell 3:  Define a function that gives the global error (see Complement 3.1, §4).

We need now to create a function that calculates the global error with reasonable efficiency. The global error is:
	e = (∑_i) { (Subscript[p, e])((Subscript[ρ, i]), T ) - pVirial[T, (Subscript[ρ, i]), B2, B3, B4] }^2. 
Here (Subscript[p, e])((Subscript[ρ, i]),T) is the measured pressure at temperature T and molar density (Subscript[ρ, i]) and pVirial[T, ρ, B2, B3, B4] is the virial equation giving the pressure at the temperature  T and molar density ρ.     I generate the global error function by using MapThread again (see Mathematica explanation below).

Clear[e]
e[f_, B2_,B3_, B4_,T_,dataList_]:= 
 Plus@@MapThread[(#1-f[T,#2,B2,B3,B4])^2&,dataList]

dataList is {pList, ρList}

Let us see if this works.

e[pVirial, B2,B3,B4,T,{pList,ρList}]

The function e[..]  does what I intended. It takes the first pressure in dataList (which is 1 atm) and subtracts from it the value calculated from the virial equation for the density 0.03498  (which corresponds to 1 atm), and squares the result.  It then adds to this, all the similar terms for all the pressures and the densities in the data lists.

Mathematica explanation:  I remind you a few Mathematica procedures. The way MapThread works is shown below.

Clear[f,A1,B1,C1,a,b,c]
	MapThread[f, {{a,b,c},{A1,B1,C1}}]

Note that f has not been defined, but because it is the first argument in MapThread the computer assumes that it is a function of two variables (because the second argument in MapThread has two lists). I can replace the function f inside MapThread with the anonymous function 
	(#1 - pVirialρ[#2, T, paramList])^2&.  
See below how this works.

MapThread[(#1-Hold[pVirial[T,#2,B2,B3,B4]])^2&,{{p1,p2,p3},{ρ1,ρ2,ρ3}}]

I used Hold[expr] to make sure that expr  is not evaluated, so that you can see what MapThread does. If I apply ReleaseHold to the above result the expression is executed

ReleaseHold[%]

Then I apply Plus to this expression to obtain the sum of the elements of the list:

Plus@@%

The one-line program defining the function e above puts into one statement  all these operations.

Cell 4:  Find the values of B2, B3, B4 that minimize the global error and design a function that gives the pressure according to the virial equation with the best coefficients we have found

I minimize the global error function with respect to B1, B2, B3. This is done in Mathematica with the function FindMinimum.  The arguments of FindMinimum are:  the function to be minimized (in our case, e) and an initial guess for the values of the virial coefficients.  The time it takes to find the minimum  is shortened by giving FindMinimum a good initial guess. I will use the poor values obtained in the WorkBook I.4. They are 
	{B2 = 0.3311, B3= -3.6081, B4 = 8.5826}.  
	
• The syntax of  FindMinimum is 
	FindMinimum[f[x,y],{x,x0},{y,y0}]
• x0 and y0 are guesses of the values of x and y where the minimum might be. 
• These guesses do not have to be accurate.

sol=FindMinimum[Evaluate[e[pVirial,B2,B3,B4,350,{pList,ρList}]],{B2,0.3311},{B3,-3.6},{B4,8.5}, MaxIterations->60]

Mathematica explanation:  FindMinimum returns the values of the coefficients B2, B3, B4 for which the error is smallest. The output contains the following information: the first element in the list, 9.33776 (10^(-7)), is the minimum value of the global error . The next element is a list of rules that give the values of the B's that produce the minimum error.  Evaluate[e[...]] forces the evaluation of the function  e  before it is used by FindMinimum.  This increases efficiency but is not necessary.

The virial coefficients are extracted by applying the rules sol[[2]] to the function pVirial

pVirial[T,v,B2,B3,B4]/.sol[[2]]

Note that the virial coefficients obtained by least square fitting are very different from the ones obtained by the crude method.

Cell 5:  Test the virial coefficients obtained by least square fitting  by comparing to experiments.

I prepare several lists that go into the table displaying the results.

The pressure calculated with the virial equation provided by the least square fitting:

pCalc=(pVirial[350,#,B2,B3,B4])& /@ρList/.sol[[2]];

In the above equation I apply the function pVirial to every density in the list ρList. Then I replace the virial coefficients with the values calculated in Cell 4. The replacement rules are named sol[[2]]

The error made in the calculation in percent:

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

Transpose the lists to put it in the format needed in TableForm:

toTableForm= Transpose[{pList, pList,pCalc,error}];

Make table:

Print["comparison with experiments at T=350K"] ;
		  TableForm[toTableForm,TableHeadings->{None,{"v L/mol","p atm","pcalc","%error"}}]

We are doing much better than in the naive fitting. The fit obtained here works better over a wider pressure range.

WorkBookT.C3.3. Determine the coefficients in the van der Waals equation

I have the following data for ethane: 
	({{"p atm", "T=300K", "T=350K", "T=400K"}, {1, 24.4300000000000006`, 28.5880000000000009`, 32.7249999999999996`}, {2, 12.1199999999999996`, 14.2289999999999983`, 16.3150000000000012`}, {3, 8.01900000000000012`, 9.44100000000000072`, 10.8450000000000002`}, {4, 5.96400000000000041`, 7.04769999999999985`, 8.11009999999999919`}, {5, 4.73299999999999965`, 5.61190000000000033`, 6.46870000000000011`}, {6, 3.9113000000000003`, 4.65470000000000005`, 5.37469999999999981`}, {7, 3.32520000000000015`, 3.97040000000000006`, 4.59330000000000016`}, {8, 2.8841000000000001`, 3.45699999999999985`, 4.00729999999999986`}, {9, 2.54210000000000002`, 3.05810000000000004`, 3.5512999999999999`}, {10, 2.26069999999999993`, 2.73870000000000013`, 3.18659999999999987`}}) 
The last three columns give the molar volumes at different temperatures and at the pressure given in column 1. 
For example, if p = 6 atm and T = 350 K, the molar volume is v = 4.6547 liter/mol.

I will use the least-squares fitting method to determine the parameters a and b in the van der Waals equation 
	p = f(T,v) = (RT/(v-b))-(a/(v^2))					(1)
so that it gives the best fit (the best that van der Waals can do)

The structure of the data is more complicated than in the previous problem.  I have to define the global error for this data.  Let me denote the pressures in the table above by (Subscript[p, i]).  The volume corresponding to the pressure (Subscript[p, i]) and the temperature (Subscript[T, α]) is denoted (Subscript[v, α,i]).  For example, if α = 2, the temperature is T = 350. If i = 3, the pressure (Subscript[p, 3]) = 3 atm.  The volume at T = 350 K and p = 3 atm is (Subscript[v, 1,3]) = 9.441 liter/mol.   For a given temperature (Subscript[T, α]) and volume (Subscript[v, α,i]) the pressure given by the van der Waals equation (see Eq. 1).is  
 	(Subscript[p, i])=((Subscript[RT, α])/((Subscript[v, α,i])-b))-(a/((Subscript[v, α,i])^2))
 The global error is (see Eq. 16 in Complement 3.1):
	e(a,b) = (∑_(α=1)^3)(∑_(i=1)^10){(Subscript[p, i])-[ ((Subscript[RT, α])/((Subscript[v, α,i]) - b))-(a/((Subscript[v, α,i])^2))]}  
The double sum appears because we sum over all the data points, three temperatures (α = 1, 2, 3) and ten volumes for each temperature ( i = 1, 2,...10). 
	To find the best values for a and b, I have to minimize  e(a, b) with respect to these two variables.

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

Cell 1:   Make lists with the data

It is convenient to organize the data in a list having the following structure. The pressures are in one list  pList = {1, 2, . . ,10}.   The volume is more complicated because for a given pressure I have different molar volumes for different temperatures. Since the data has only three temperatures, I make three lists named v[1], v[2], v[3] which contain the molar volumes at the three temperatures T = 300K, T = 350 K and T = 400 K.
	Here is how you make contact with the notation used earlier:  The index α labels v[α]; the index i gives the position of the volume in the list v[α]. For example, if α = 2 and i = 6, then v[2][[6]] is the volume (Subscript[v, α,i ])= 4.6547 liter/mol.

pList=Table[i,{i,1,10}] (*pressures in atm*);

At T = 300 K the molar volumes at the pressures in pList are (see Appendix 2):

v[1]  (*molar volume at 300K*)=
	{24.430,12.12,8.019,5.964,4.733,3.9113,3.3252,2.8841,2.5421,2.2607}(*liter/mol*) ;

At T = 350 K the molar volumes at the pressures in pList are (see Appendix 2):

v[2](*molar volume at 350 K*)={28.588,14.229,9.441,7.0477,5.6119,4.6547,3.9704,3.457,3.0581,2.7387}(*liter/mol*);

At T = 400 K the molar volumes at the pressures in pList are (see Appendix 2):

v[3](*molar volume at 400K*)={32.725,16.315,10.845,8.1101,6.4687,5.3747,4.5933,4.0073,3.5513,3.1866} (*liter/mol*);

Cell 2:  A function that calculates the pressure with the van der Waals equation.

Define pvdW[v,T] which calculates the pressure from the van der Waals equation p=((R T)/(v - b ))-(a/(v^2)).

Clear[a,b,pwdW];
pvdW[T_,v_]:=With[{R=0.082067 (*liter atm/K mol*)},
Clear[a,b];
		(R T/(v-b)-a/v^2)];
SetAttributes[pvdW,HoldAll]

Here is how this function is used:

pvdW[T,v ]

Cell 3  Define a procedure that calculates the error e(a,b).

Calculate the errors for each temperatures and then add them all up.

e[a_,b_]:=Module[{s1,s2,s3},
s1=MapThread[(#1-pvdW[300,#2] )^2&,{pList,v[1]}]; 
(*s1 is a list of the errors squared, for T=300K*)
s2=MapThread[(#1-pvdW[350,#2])^2&,{pList,v[2]}];
(*s2 is a list with the error squared, for T=350 K*) 
  s3=MapThread[(#1-pvdW[400,#2])^2&,{pList,v[3]}]; 
(*s3 is a list of errors squared, for T=400K*)
(*below: use Join to make a list with all square errors and then apply Plus@@  to add them*) 
Plus@@(Join[s1,s2,s3])]

Use the function Short to look at a few terms in e[a,b].    
	e[a,b]//Short[#,7]& 
controls the output and I get only a few of the terms in the expression.  I can see that the function works as I intended

e[a,b]//Short[#,7]&

Cell 4  Find the values of a and b that minimize the global error.

sol=FindMinimum[Evaluate[e[a,b]],{a,0.027},{b,1.358},
	MaxIterations->60]

The second term in the list sol is the list of two rules that will replace a and b with the numerical values that minimize the errors.  Here is how you calculate the error

e[a,b]/.sol[[2]]

The expression below returns the van der Waals formula for  pressure with numerical values for a and b that I determine to minimize the error.

pvdW[T,v]/.sol[[2]]

Cell 5:  Test how well the result fits data.

Below I calculate three lists vdWPressure[i], i = 1, 2, 3, containing the pressure given by the van der Waals equation for the three temperatures (e.g. vdWPressure[1] gives the pressures for T=300 K and the molar volumes in the list vList[1]).

vdWPressure[1]=(pvdW[300,#] &/@v[1])/.sol[[2]];
vdWPressure[2]=(pvdW[350,#]&/@v[2])/.sol[[2]];
vdWPressure[3]=(pvdW[400,#]&/@v[3])/.sol[[2]];

Organize the elements of the lists so that you can print them with TableForm:

tb1=MapThread[{300,#1,#2,#3,(#2-#3)*100/#2}&,{v[1],pList,vdWPressure[1]}];(Print["comparison with experiments at T = 300K/n"];TableForm[tb1,TableHeadings->{None,{"T","(Subscript[v, exp])","(Subscript[p, exp])","(Subscript[p, vdw])","%error"}}])

tb2=MapThread[{350,#1,#2,#3,(#2-#3)*100/#2}&,{v[2],pList,vdWPressure[2]}];(Print["Comparison to experiment at T = 350K "]; TableForm[tb2,TableHeadings->{None,{"T K","(Subscript[v, exp])","(Subscript[p, exp])","(Subscript[p, vdw])","%error"}}])

tb3=MapThread[{400,#1,#2,#3,(#2-#3)*100/#2}&,{v[3],pList,vdWPressure[3]}];(Print["Comparison to experiment at T = 400K/n"];TableForm[tb3,TableHeadings->{None,{"T K","(Subscript[v, exp])","(Subscript[p, exp])","(Subscript[p, vdw])","%error"}}])

Cell 6:  Compare to the van der Waals equation given in textbooks.

Appendix 4.1 gives for ethane the coefficients shown below (I converted bar to atm):

coeffsTextBook={a->5.4897,b->0.06297};

I calculate the pressure given by the van der Waals equation with the textbook coefficients to see if it does better or worse than mine when compared to experiments. Use the data for 400 K.

pressureVdWTextB400=(pvdW[400,#]/.coeffsTextBook)&/@v[3]; Print["The temperature T = 400 K"] ;
	TableForm[MapThread[{#1,#2,#3,(#2-#3)*100/#2,#4,(#2-#4)*100/#2}&,{v[3],pList,vdWPressure[3],pressureVdWTextB400}],
	TableHeadings->{None,{"v","p","pmine","errorMine","pText","ErrorTex"}} ]

Look how the equations work at T = 300K

pressureVdWTextB300=(pvdW[300,#]/.coeffsTextBook)&/@v[1]; Print["The temperature T = 400 K"] ;
	TableForm[MapThread[{#1,#2,#3,(#2-#3)*100/#2,#4,(#2-#4)*100/#2}&,{v[1],pList,vdWPressure[1],pressureVdWTextB300}],
	TableHeadings->{None,{"v","p","pmine","errorMine","pText","ErrorTex"}} ]

It is clear that the coefficients I determined work better than the ones I found in the literature. This happens because I fitted the data in the range in which I tested the equation.