WorkBook T11. The Change of Entropy

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

WorkBookT. 11.1  The isobaric change of entropy

Cell 1: Define the problem, give the data

• Calculate the change of entropy of a mole of (Subscript[N, 2]),  in the transformation 
	(Subscript[σ, i]) = { (Subscript[p, i]) = 1 atm, (Subscript[T, i]) = 300 K} -> (Subscript[σ, f]) = { (Subscript[p, f]) = 1 atm, (Subscript[T, f]) = 800 K}
• From literature I have (for p = 1 atm)
	(Subscript[c, p]) = 26.924 + 5.910 ×(10^(-3)) T - 3.376 ×(10^(-7)) (T^2) cal/mol K	  (1)
 • Compare the result with the tabulated experimental entropies. At the temperatures
 	 T = {300, 400, 500, 600, 700} K 
 the entropies of (Subscript[N, 2]) are
 	s = {191.37, 199.72, 206.22, 211.56, 216.08} cal/mol K.
 •  (Subscript[c, p]) and the entropies are give in cal/mole K and will work with these units. 
 • If agreement between the calculated and the tabulated entropy is poor, determine the cause of disagreement and fix the problem.

Input the data:

sDat ={191.37,199.72,206.22,211.56,216.08}(*cal/mol K*)

Define a function that gives
	(Subscript[c, p]) =  26.924 + 5.910 ×(10^(-3)) T - 3.376 ×(10^(-7)) (T^2)     J/mol K	 (1)

Clear[T,cp]	
cp[T_]:=26.924 + 5.910 (10^-3) T -3.376 (10^-7 )T^2 (*cal/Kmol K*)

cp[500]

Cell 2: Calculate the change of entropy

I have at my disposal two equations:
	ds = ((Subscript[c, v])/T) dT +  (Subscript[(∂p /∂T), v]) dv         					(2)
or 
	ds = ((Subscript[c, p])/T) dT - (Subscript[(∂v/∂T), p]) dp 						(3)
For an isobaric transformation dp = 0 so Eq. 3 is simplest. Therefore, the variation of the entropy with temperature, when the pressure is held constant, is:
Δs(T) ≡ s(T,p=1)-s(300,p=1) = (∫_( 300)^T) (((Subscript[c, p])(T, p=1))/T)dT 	  	(4) 
This is Eq. 11.17 in the textbook.  It gives the change in entropy in cal/mole K, for a transformation in which the temperature changes from 300 K to  T  K. The units are cal/mole K because (Subscript[c, p]) is in cal/mole K

Clear[Δs]
Δs[T_]= Integrate[cp[T]/T,{T,300,T},GenerateConditions->False] (* J/mol K*)

The entropy change ought to be of order  tens of Joule

Δs[600]

Entropy must increase with temperature

Δs[600]>Δs[500]> Δs[300]

Cell 3: Compare the calculated Δs to the values given in literature

I will compare the results of my calculation to the entropy values given in the literature.  I found data that says that at p = 1 atm and 
	T = {300, 400, 500, 600, 700} K 
 the entropy of (Subscript[N, 2]) is
 	s = {191.37, 199.72, 206.22, 211.56, 216.08} cal/mol K.
I will now check whether my calculations give results close to these.

• There is a slight complication that I must take care of. 
• I can calculate an entropy difference
	 s(T, p = 1) - s(T = 300, p = 1) = (∫_( 300)^( T )) ((Subscript[c, p(T,p=1)])/T)dT		(1)
for any value of T.   
• However, the table with data gives me the entropy  (not an entropy difference).  
• So I need to use the entropies from the table to calculate the entropy differences 
	{s(400) - s(300), s(500) - s(300), s(600) - s(300), s(700) - s(300)}
• I can compare these differences, calculated by using the entropy from the table, to the differences calculated with Eq. 1.

The entropy differences calculated with the data from the tables are

sDifDat=sDat-sDat[[1]] (*Joule/mole K*)

I need to get rid of the first term in the list

sDifDat=Drop[ sDifDat,1]

Now sDifDat is  
	{s(400) - s(300), s(500) - s(300), s(600) - s(300), s(700) - s(300)}

The temperatures for which I have entropy differences are

Temp=Table[T,{T,400,700,100}]

Now calculate the   difference 
	s(T, p = 1) - s(T = 300,p = 1) =  (∫_300^T)(Subscript[c, p])(T,p=1 atm)/TdT  
with the heat capacity given by Eq. (2) for all the temperatures given in the list Temp

ours =  Integrate[cp[T]/T,{T,300,#}]&/@Temp(*cal/mol K*)

Mathematica explanation: the function Integrate[cp[T]/T, {T, 300 ,#}]& is mapped (symbol /@) into the list of temperatures TempList. This calculates a list of  integrals:  the final temperature in each integral is one of the temperatures in the list TempList.

Make a table with the results of interest

TableForm[
	MapThread[{#1,#2,#3,100*(#2-#3)/#2}&,{Temp,sDifDat,ours}],
	TableHeadings->{None,{"T K","ΔSdat","ΔScalc","%error"}}]

The error is not acceptable. I need to find the reason for the discrepancy.

Cell 4:  Test whether the discrepancy is due errors in (Subscript[c, p])

• I suspect that the error comes from the equation used for (Subscript[c, p]) (which I got from literature).  I will create my own function (Subscript[c, p])(T) by fitting a large set of good quality data.
• For 
	p = 1 atm,  
and for
 	T = {300, 320, 340, 360, 380, 400, 450, 500, 550, 600, 650, 700} 
Din’s book gives (for (Subscript[N, 2])):
	(Subscript[c, p]) = { 28.98, 29.00, 29.02, 29.05, 29.07, 29.09, 29.15, 29.21, 29.27, 29.33, 29.39, 29.45} joules/mol K. 
• Let us check if values of (Subscript[c, p]) calculated with equation (2) are close to the values given by Din.

Make a list with the temperatures at which (Subscript[c, p]) data is available and a list with the values of (Subscript[c, p]) (from Din) at those temperatures.
I put the temperatures for which Din gives (Subscript[C, p]) in the list TempCp and the experimental values of (Subscript[C, p]) in the list CpDat

TempCp={300,320,340,360,380,400,450,500,550,600,650,700}; 
CpDat={ 28.98,29.00,29.02,29.05,29.07,29.09,29.15,29.21,29.27,29.33,29.39,29.45}(*Joule/K mole*);

Calculate (Subscript[c, p]) given by Eq. (2) (which I suspect to have errors) for the temperatures in the list TempCp for which I have data:

CpCalc=  cp[#]&/@TempCp

Make a table with these values and the error made by my equation for (Subscript[c, p])

TableForm[
	MapThread[{#1,#2,#3,100*(#2-#3)/#2}&,{TempCp,CpDat,CpCalc}],
TableHeadings->{None,{"T K","Cpexp","CpCalc","%error"}}]

• The equation for (Subscript[c, p]) that I used (which is from the literature: a source other than Din) does not fit the data very well.  
• This can  be the reason for the difference between the calculated entropy and that given by Din. 
• I create a function (Subscript[c, p])(T)  that fits Din's data for (Subscript[c, p]) and then use it to calculate the change in entropy with temperature.  
• If my suspicion is correct this should give better values for the entropy changes.

I use Mathematica function Fit to fit the data for (Subscript[c, p]) (from Din) to a polynomial of third order in T.

Clear[CpFit]
CpFit[T_]=Fit[MapThread[{#1,#2}&,{TempCp,CpDat}],{1,T,T^2,T^3},T]  (*cal/mol K *)

Mathematica explanation:  I want to fit the data with a polynomial of the form 
	(Subscript[c, p]) = a + b T + c (T^2) + d (T^3). 
This means that I want to get the values of a, b, c, for which the polynomial reproduces best the data stored in CpDat. Mathematica provides for this a function called Fit. The input to Fit is a list of the form 
	{{var1,data1},{var2,data2},...} 
where var1 is the first value of temperature and data1 is the corresponding value of (Subscript[c, p]), etc. I put my data in this form with 
	MapThread[{#1, #2}&,{TempCp,CpDat}]

a=MapThread[{#1,#2}&,{TempCp,CpDat}]

Then I give  this information to Fit:

Fit[a,{1,T,T^2,T^3},T]

Test the function

CpFit[T]

Compare to the function given by Eq. (2)

cp[T]

The two expressions are  different.  Check the goodness of the fit by comparing the function with the data fitted

TableForm[MapThread[{#1,#2,#3,(#2-#3)*100/#2}&,{TempCp,CpDat,CpFit/@TempCp}],TableHeadings->{None,{"T K","Cpexp","CpFit","fitErr"}},
		TableSpacing->{1,2}]

As you can see my fit of the data is better than that provided by Eq. (2).  Let us use the CpFit to calculate the change of entropy and compare the result to Din's.

Cell 5:  Calculate the change of entropy with Eq. (1), but use the function for (Subscript[c, p]) determined by me by fitting Din's (Subscript[c, p]) data

The calculated entropy change is (this is s(T) - s (300) ):

Temp

SfromFit= Integrate[CpFit[T]/T,{T,300,#}]&/@Temp

TableForm[MapThread[{#1,#2,#3,100*(#2-#3)/#2}&,
			{Temp,sDifDat,SfromFit}],
	TableHeadings->{None,{"T K","ΔSdat","ΔScalc","%error"}}]

The agreement between the entropy calculated from the (Subscript[c, p]) obtained by fitting Din's (Subscript[c, p]) data and the entropies tabulated by Din is very good. 

You learn from this that you must check any equation obtained from the literature, if the calculation is important. Often there are errors in the literature (typos or mistakes in measurements or in calculations) and you should not assume that everything that has been printed is correct.

WorkBookT. 11.2 The isothermal change of entropy

Cell 1:  Define the problem

• Calculate the change of entropy when 3 moles of (Subscript[CO, 2]) are compressed from (Subscript[p, i]) =1 atm to (Subscript[p, f]) =100  atm at a constant temperature equal to (Subscript[T, f]) = 1000 K. 
• Use the Benedict-Webb-Rubin equation of state. 
• Then use the ideal gas law and compare the entropy calculated this way with the one obtained by using the BWR equation.

Use the equation 
 s(T = 1000K, p = 100 atm) - s(T = 1000K, p = 1atm) = 	(∫_( (Subscript[v, i]))^( (Subscript[v, f]))) (Subscript[(((∂p)/(∂T))), v])dv  	
The difficult part is calculating the initial and final molar volumes.

Cell 2:  Define a function that gives the Benedict Web Rubin equation of state

Use the Benedict Webb Rubin equation defined in a previous WorkBook

Clear[pBWR]
pBWR[T_,v_,parameters_]:=Module[
{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];

Check that the equation is correct

pBWR[T,v,{}]

Get the parameters in BWR equation for (Subscript[CO, 2]) (from Appendix 5)

param={A->2.7374, B->4.9909*10^-2,C0->1.38564*10^5,a->1.3681*10^-1,b->7.2105*10^-3,c->1.49183*10^4,α->8.4658*10^-5,γ->5.393*10^-3,R->0.08205};

As a check the function PBWR I perform the following calculation: at  T = 500 K  and a molar volume of  40 liter/mol the gas should be ideal. Check if the function PBWR gives the same result as the ideal gas law

pBWR[400,40,param]

Use ideal gas law to check if this number is reasonable.

0.082 400/40

They are close enough to  convince me that I use the function PBWR correctly and the function works  well

Cell 3:  The derivative (Subscript[( ((∂p)/(∂T))), v])

D[pBWR[T,v,{}],T]//Expand

Cell 4: The integral  ∫(Subscript[(((∂p)/(∂T))), v])dv

Integrate[D[pBWR[T,v,{}],T],v]//Expand

Cell 5:  Calculate the initial and final molar volumes

I calculate the molar volumes by using the function VolumeG defined in WorkBook3.   The equation of state is Benedict-Webb-Rubin.

Copy below the function volumeG:

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}] ];

Use volumeG to obtain replacement rules fot the initial and final volumes

volRules={vi->volumeG[1000,1,pBWR,param],vf->volumeG[1000,100,pBWR,param]}

The initial volume is 82.073 liter/ mol and the final volume is 0.842 liter/mol.

Cell 6:  Calculate the change of entropy

Use
	 Δs = ∫(Subscript[(((∂p)/(∂v))), T])dv  
The equation gives the change of entropy  in liter atm/ K mol. Convert it to joules by using
	1 liter atm = 101.325 J

Calculate  
	Δs =(∫_( (Subscript[v, i]))^(Subscript[v, f]))(Subscript[(((∂p)/(∂T))), v])dv    for constant T

Δs=101.325*Integrate[D[pBWR[T,v,param],T]/.T->1000, {v,vi,vf},GenerateConditions->False]/.volRules     (*gives the entropy change in J/mol K*)

For three moles the change in entropy is

3*Δs (*J/K *)

Cell 6:  Is the error made by using the ideal gas law large?

Calculate the entropy for the ideal gas from:
	Δs = R ln  {((Subscript[v, f])/(Subscript[v, i]))} = R ln{((Subscript[p, i])/(Subscript[p, f]))}

ΔsIdeal=3* 8.31451 * Log[1/100]

This is amazing.  The calculation using BWR gave -115.912 J/K.  Even though we go to 100 atm the ideal gas law works fairly well as far as the change in entropy is concerned.  This happens because the temperature is high and the ideal gas law is OK. In the next cell I show that at lower temperature the ideal gas law gives a big error in the entropy change.

Cell 7:  Entropy change when p  changes from 1 atm to 100 atm and T = 298 K.

I calculate the isothermal transformation at T = 298 K when the ideal gas law should not work.

Calculate the molar volumes given by the BWR equation, at 298 K are and p = 1 atm or p =100 atm.

volRules={vi->volumeG[298,1,pBWR,param],vf->volumeG[298,100,pBWR,param]}

Calculate the entropy from  
	Δs =(∫_( (Subscript[v, i]))^(Subscript[v, f]))(Subscript[(((∂p)/(∂T))), v])dv 
 The equation gives the entropy in liter atm/ K mol. Convert it to joules by using
	1 liter atm = 101.325 J

Δs=101.325*Integrate[D[pBWR[temperature,v,param],temperature]/.temperature->298, {v,vi,vf}]/.volRules     (*gives the entropy change in J/mol K*)

Calculate the entropy per mole, for the same transformation, by using the ideal gas law.  The formula for entropy is (you obtain this by integrating Eq. 2 and using the ideal gas law to obtain the integrand)
	Δs = R ln  {((Subscript[v, f])/(Subscript[v, i]))} = R ln{((Subscript[p, i])/(Subscript[p, f]))}

ΔsIdeal= 8.31451 * Log[1/100] (*J/K mol*)

You see that at 298 K the result obtained by using the ideal gas law is very different from that obtained by using the BWR equation.

WorkBookT. 11.3. The independence of entropy on transformation paths.

Cell 1:  Define the problem

Calculate the entropy change for a mole of (Subscript[CH, 4]) in the transformation  
	(Subscript[σ, i])={(Subscript[T, i])=298 K, (Subscript[p, i])=1 atm}  -> 	
	(Subscript[σ, f])={(Subscript[T, i])=1000  K, (Subscript[p, i])=100 atm}, 
I will perform the calculation along two paths to verify that the entropy change is independent on the path (all paths lead to the same entropy change  as long as initial and the final states are the same).

Path I consists of two segments 
Path  IA:  (Subscript[σ, i])={ 298 K, 1 atm} -> (Subscript[σ, 1])= {298K, 100 atm}
Path IB:  (Subscript[σ, 1])= {298K, 100 atm}-> (Subscript[σ, f]) = {1000 K, 100 atm}
The entropy change for (Subscript[σ, i]) -> (Subscript[σ, f]) is
	 (Subscript[Δs, IA]) + (Subscript[Δs, IB])

Path II consists of two segments 
Path  IIA:  (Subscript[σ, i])={ 298 K, 1 atm} -> (Subscript[σ, 2])= {1000K, 1 atm}
Path IIB:  (Subscript[σ, 2])= {1000K, 1 atm}-> (Subscript[σ, f]) = {1000 K, 100 atm}
The entropy change for (Subscript[σ, i]) -> (Subscript[σ, f]) along path B is
	(Subscript[Δs, II]) = (Subscript[Δs, IIA]) + (Subscript[Δs, IIB])

Theory says that entropy is path independent; it depends only on the initial and the final state.  If this is the case I will find that 
	(Subscript[Δs, IA]) + (Subscript[Δs, IB]) = (Subscript[Δs, IIA]) + (Subscript[Δs, IIB])
We encounter here a very tricky situation.  Second law says that without doubt the change of entropy does not depend on path. We want to verify this by numerical calculations.  For this we must calculate the entropies with perfect accuracy. Any approximation ruins the equality.

Cell 2: Collect the necessary information

I need the heat capacity of methane and a good equation of state. 
I use (from Appendix 8):
(Subscript[C, p])(T, p=1) = 4.171 + 14.450×(10^(-3))T + 1-0.267×(10^(-6))(T^2) cal/mol K	   (2)
If T is in K the result for (Subscript[c, p]) is in cal/mol K

Define a function that calculates (Subscript[c, p])

Clear[cp]
cp[T_]= 4.171 + 14.450 10^-3 T -0.267 10^-6 T^2 (*cal/ mole K*)

Check that this gives the right order of magnitude:

cp[300](*cal/mole K*)

The change of entropy with the temperature when the pressure is held constant and equal to (Subscript[p, 0]) is given by
	(Subscript[(Subscript[Δs, p]),  ]) = s((Subscript[T, f]), (Subscript[p, 0])) - s((Subscript[T, i]),(Subscript[p, 0])) = 4.184 (∫_( (Subscript[T, i]))^(Subscript[T, f]))(Subscript[c, p])(T,(Subscript[p, 0])) ddT       Joule
I don’t have the dependence of (Subscript[c, p]) on pressure and I assume that (Subscript[c, p]) is pressure independent. This is a bad approximation as far as methane is concerned.

The change of entropy when the volume changes and the temperature is held constant is
	(Subscript[Δs, T])   =  s(T,(Subscript[v, f])) - s(T,(Subscript[v, i])) = 101.325 (∫_( (Subscript[v, i]))^( (Subscript[v, f]))) (Subscript[(((∂ p(T,v))/(∂T))), v])dv        Joule

To calculate ∂p(T,v)/∂T,  (Subscript[v, i]) and (Subscript[v, f]) I need an equation of state.  I will use the Benedict, Webb and Rubin (BWR) equation (see Chapter 2) which is
	p = R T ρ + (B R T- A-(Subscript[C, 0])/(T^2))(ρ^2)+ (b R T -a) (ρ^3)+a α (ρ^6)+
	      c (ρ^3)(1 + γ (ρ^2)) exp[-γ (ρ^2)]/(T^2)

Define the Benedict-Webb-Rubin equation

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];

Define the replacement rules for methane.  For (Subscript[CH, 4]) the parameters are:
A = 1.855,  B = 0.0426,  C0 = 0.02257 × (10^6), a = 0.0494,  b = 0.00338,  c = 0.002545 (((×10))^6), γ = 0.006,  α = 0.00012436 
(From Din, Thermodynamic Functions for Gases) The units are such that inserting ρ =1/v in mole/liter, T in K, and R = 0.0820 liter atm/mol K gives the pressure in atmospheres.

parametersBWR={A->1.855,B->0.0426, C0->0.02257 10^6, a -> 0.0494,b->0.00338, c->0.002545 10^6,γ->0.006,α->0.00012436}

Test that the pressure is close to 1 atm when T = 298 K and v = 24 liter/mole

pBWR[298,24,parametersBWR]

Cell 3: Calculate the change of entropy along path I.

Path I consists of two segments 
Path  IA:  (Subscript[σ, i])={ 298 K, 1 atm,(Subscript[v, i])} -> (Subscript[σ, 1])= {298K, 100 atm,(Subscript[v, 1])}
Path IB:  (Subscript[σ, 1])= {298K, 100 atm,(Subscript[v, 1])}-> (Subscript[σ, f]) = {1000 K, 100 atm,(Subscript[v, f])}
The entropy change when going from (Subscript[σ, i]) to (Subscript[σ, f]) is along path I is 
(Subscript[Δs, IA]) + (Subscript[Δs, IB]) =   101.325 (∫_( (Subscript[v, i]))^(Subscript[v, 1]))(((∂p(T, v))/(∂T)))(Subscript[|, T=298])dv  + 4.184 (∫_( (Subscript[T, i]))^( (Subscript[T, f])  ))((Subscript[c, p(T, p=100)])/T)dT 
I want the units to be joule/K mole and therefore I use 
	 1 liter atm  = 101.325 Joule
and 
	1 cal = 4.184 Joule
To evaluate the integral over the volume I need to calculate  (Subscript[v, i])  and (Subscript[v, 1]) by solving 
	1atm = f(298 K, (Subscript[v, i]))
	1 atm =  f( 1000,(Subscript[v, 1]))   
	100 atm = f (298, (Subscript[v, f]))
I do this by using the function volumeG defined earlier and reproduced below

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}] ];

volumeRules={vi->volumeG[298,1,pBWR ,parametersBWR],
v1->volumeG[298,100,pBWR ,parametersBWR],vf->volumeG[1000,100,pBWR ,parametersBWR] }

Evaluate the entropy along the path I:
(Subscript[Δs, IA]) =  101.325 (∫_( (Subscript[v, i]))^(Subscript[v, 1])) (((∂p(T, v))/(∂T))) dv /.{T->298 ,v1->0.2075,(Subscript[v, i])->24.4073}
(Subscript[Δs, IB]) =  4.184 (∫_( 298)^( 1000  ))((Subscript[c, p(T, p=100)])/T)dT

Here I made an error. I need (Subscript[c, p])(T, p = 100)  but I don’t know it and I use the  
(Subscript[c, p])(T, p = 1) instead.

ΔsIA = 101.325*Integrate[D[pBWR[T,v,parametersBWR],T],{v,vi,v1},GenerateConditions->False]/.volumeRules/.T->298

Increasing the pressure decreases the entropy always, because
	 (Subscript[(((∂s)/(∂p))), T]) = - (Subscript[(((∂v)/(∂T))), p])
and the right-hand-side is positive (the volume increases when I increase the pressure)

The entropy change on path IB is

ΔsIB=4.184Integrate[cp[T]/T,{T,298,1000}](*Joules/mole K*)

Increasing the temperature, at constant pressure, always increases the entropy because the derivative of the entropy with temperature is positive (it is (Subscript[c, p])/T).

The entropy change on path I is

ΔsI = ΔsIA+ΔsIB   (*Joule/K mole*)

Cell 4: Calculate the change of entropy along path II.

Path II consists of two segments 
Path  IIA:  (Subscript[σ, i])={ 298 K,1 atm,(Subscript[v, i])} -> (Subscript[σ, 2]) = {1000K, 1 atm, (Subscript[v, 2])}
Path IIB:  (Subscript[σ, 2]) = {1000K, 1 atm,(Subscript[v, 2])} -> (Subscript[σ, f]) = {1000 K, 100 atm, (Subscript[v, f])}
The entropy change for the transformation (Subscript[σ, i]) -> (Subscript[σ, f]) along path B is
	(Subscript[Δs, II]) = (Subscript[Δs, IIA]) + (Subscript[Δs, IIB])

with
	(Subscript[Δs, IIA]) = 4.184 (∫_( 298)^( 1000)) (((Subscript[c, p])(T, p=1))/T)dT       Joule
and
	(Subscript[Δs, IIB]) = 101.325 (∫_( (Subscript[v, 2]))^(Subscript[v, f]))(((∂p(T, v))/(∂T)))(Subscript[|, T=1000])dv      Joule

ΔsIIA = 4.184 Integrate[cp[T]/T,{T, 298,1000}] (*Joule*)

To calculate (Subscript[Δs, IIB]) I need to calculate (Subscript[v, i]), (Subscript[v, 2]) and (Subscript[v, f]) by solving
	1 = f(298,(Subscript[v, i]))
	1 = f(1000,(Subscript[v, 2]))
	100 = f(1000,(Subscript[v, f]))

Clear[volumeRules]
volumeRules={vi->volumeG[298,1,pBWR ,parametersBWR],v2->volumeG[1000,1,pBWR ,parametersBWR],
vf->volumeG[1000,100,pBWR ,parametersBWR]}

ΔsIIB=101.325 Integrate[D[pBWR[T,v,parametersBWR],T],{v,vi,vf},GenerateConditions->False]/. volumeRules /.T->1000

ΔsII =ΔsIIA+ΔsIIB

{ΔsI,ΔsII}

Theory says that these to quantities must be equal and they are not.  I expected that.  If I want to test that two quantities are equal.  I can’t do that if I compute them approximately.  What we discovered is that we made some approximations; if we didn’t these two quantities would be equal.

WorkBookT. 11.4.  The independence of entropy on path: Ideal gas

Cell 1.  Collect data and equations

Path I consists of two segments 
Path  IA: 
	 (Subscript[σ, i])={ 298 K, 1 atm,(Subscript[v, i])} -> (Subscript[σ, 1])= {298K, 100 atm,(Subscript[v, 1])}
Path IB: 
	 (Subscript[σ, 1])= {298K, 100 atm,(Subscript[v, 1])}-> (Subscript[σ, f]) = {1000 K, 100 atm,(Subscript[v, f])}
The entropy change when going from (Subscript[σ, i]) to (Subscript[σ, f]) is along path I is 
	(Subscript[Δs, IA])  =  (∫_( (Subscript[v, i]))^(Subscript[v, 1]))(((∂p(T, v))/(∂T)))dv /.T->298  
	(Subscript[Δs, IB])  = (∫_( (Subscript[T, i]))^( (Subscript[T, f])  ))(((Subscript[c, p])(T, p=100))/T)dT 
	(Subscript[Δs, I]) = (Subscript[Δs, IA]) + (Subscript[Δs, IB])
I need to pay attention to units when I perform evaluations.

Define the equation of state and (Subscript[c, p])

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

Clear[pIdeal]
pIdeal[T_,v_]:=R T/v
pIdeal[T,v]
cp[T_]=5 R/2

The entropy change on the isothermal path IA
	(Subscript[Δs, IA])  =  (∫_( (Subscript[v, i]))^( (Subscript[v, 1])))(((∂p(T, v))/(∂T)))dv /.T->298

Cell 2.  Entropy change on path I

(Subscript[Δs, IA])  =  (∫_( (Subscript[v, i]))^(Subscript[v, 1]))(((∂p(T, v))/(∂T)))dv /.T->298

ΔsIA= Integrate[D[pIdeal[T,v],T],{v,vi,v1},GenerateConditions->False]

In this equation R = 8.3145 Joule.  Log[(Subscript[v, 1])/(Subscript[v, i])] is dimensionless

Off[Solve::ratnz]

Calculate the volume change on path I

volumeRules= { Solve[1==pIdeal[298,vi],vi] ,Solve[100==pIdeal[298,v1]],Solve[100==pIdeal[1000,vf]]} /.R->0.082//Flatten

ΔsIA=ΔsIA/.R->8.3145/.volumeRules

Calculate
	(Subscript[Δs, IB])  = (∫_( (Subscript[T, i]))^( (Subscript[T, f])  ))(((Subscript[c, p])(T, p = 100))/T)dT

ΔsIB=Integrate[(5 R/2)/T,{T,Ti,Tf},GenerateConditions->False]/.Tf->1000/.Ti->298/.R->8.3145

The entropy change on path I is:

ΔsI= ΔsIA +ΔsIB

Cell 3.  Entropy change on path II

Look at path II. 
path IIA:	{298, 1, (Subscript[v, i])} -> {1000, 1, (Subscript[v, 2])} 
path IIB: {1000,1,(Subscript[v, 2])} ->{1000, 100,(Subscript[v, f])}

(Subscript[Δs, IIA]) = (∫_( 298)^( 1000)) (((Subscript[c, p])(T, p=1))/T)dT      
and
	(Subscript[Δs, IIB]) = (∫_( (Subscript[v, i]))^(Subscript[v, f]))(((∂p(T, v))/(∂T)))(Subscript[|, T=1000])dv
The units will be controlled by the units of R.

Calculate 
	(Subscript[Δs, IIA]) = (∫_( 298)^( 1000)) (((Subscript[c, p])(T, p=1))/T)dT

ΔsIIA= Integrate[(5 R/2)/T,{T,298,1000},GenerateConditions->False]

ΔsIIA =ΔsIIA/.R->8.3145    (*Joules*)

Calculate
	(Subscript[Δs, IIB]) = (∫_( (Subscript[v, i]))^(Subscript[v, f]))(((∂p(T, v))/(∂T)))(Subscript[|, T=1000])dv

ΔsIIB=Integrate[D[pIdeal[T,v],T],{v,v2,vf},GenerateConditions->False]

volumeRules= { Solve[1==pIdeal[298,vi],vi] ,Solve[1==pIdeal[1000,v2]],Solve[100==pIdeal[1000,vf]]} /.R->0.082//Flatten

ΔsIIB = ΔsIIB/.{R->8.3145}/.volumeRules (*Joules*)

ΔsII=ΔsIIA +ΔsIIB

Cell 4. Compare entropy on path I to entropy on path II

{ΔsI,ΔsII}

WorkBookT. 11.4.  Optional:  A study of an adiabatic transformation

Cell 1:  Calculate how the temperature changes with the  molar volume in an adiabatic transformation.

The gas is acetylene, the initial temperature is T=300 K, and the initial volume is (Subscript[v, i]) = 22 liter/mol. Assume ideal gas and (Subscript[c, v]) independent of temperature.

Use the equation
	T = (((Subscript[T, i])((Subscript[v, i])/v))^(R/(Subscript[c, v])))		(1)

Define a function that evaluates the temperature according to Eq. 1. (Subscript[c, v]) for acetylene = 8.584 J/mol K. Remember that in deriving this equation in Supplement 11.3, I  assumed that the gas is ideal and that (Subscript[c, v]) is independent of temperature and pressure.

See below a function that calculates the temperature at the end of an adiabatic transformation performed on acetylene, when the initial temperature is 300 kelvin and the initial volume is 22 liter/mol

tempAdiab[v_]:= Module[{R=1.987 (*cal/mol K*),vi=22 (*liter/mol*),
			Ti=300 (*K*),cv=8.584 (*cal/mol K*)},
		T=Ti*(vi/v)^(R/cv)]

tempAdiab[v]

Make a plot.

Plot[tempAdiab[v],{v,22,0.3},PlotRange->All,
AxesLabel->{"v liter/mol","T K"}]

Cell 2: Calculate the dependence of pressure on the molar volume during an adiabatic transformation.

Define a function that give pressure as a function of volume and temperature, for an ideal gas:

p[v_,T_]:=With [{R=0.082 (*liter atm/mol K*)},
			R T/v]

In an adiabatic transformation the temperature is a function of the volume. The function tempAdiab[v], defined in Cell 1, gives this dependence.  Use this when you calculate the pressure along the adiabatic transformation:

p[v,tempAdiab[v]]

Plot[p[v,tempAdiab[v]],{v,10,0.3}, PlotRange->All,
AxesLabel->{"v liter","p atm"}]

I can make a plot of the temperature versus pressure in an adiabatic transformation

tb=Table[{ p[v,tempAdiab[v]],tempAdiab[v]},{v,0.2,22,0.2}] ;

Clear [T,p]
ListPlot[tb, AxesLabel->{"p atm","T K"},
PlotRange->All, PlotStyle->PointSize[0.02]]

At pressures of over 350 atm, the temperature approaches 900 K.  Keep in mind that we have made some very crude assumptions (ideal gas law, (Subscript[c, v]) independent of temperature and pressure).  If you feel brave, try removing these assumptions and see how much more complicated the problem becomes. Nevertheless it is manageable if you stay cool and know your math (you will have to solve a differential equation by using NDSolve).