© 2022 Horia Metiu

WorkBookT.14:  The change of chemical potential in a thermodynamic transformation

WorkBookT.14.1. The change of chemical potential in an isobaric transformation

Calculate the change in chemical potential in the transformation
{(Subscript[T, i])=298.15 K, p=1 atm} ⟶ {(Subscript[T, f])=130 K, p =1 atm}

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

Plan: I will use the equation (14.6 in textbook)
	μ((Subscript[T, f]),(Subscript[p, i])) -μ((Subscript[T, i]),(Subscript[p, i])) = - (∫_(Subscript[T, i])^(Subscript[T, f]))s(T, (Subscript[p, i]))dT		 (1)
To evaluate the integral I need to calculate how the entropy 
	s(T, p = 1 atm) 
depends on temperature. I use Eq. 14.9 from the textbook:
s(T, (Subscript[p, i]) ) =  s(T = 298.15 , (Subscript[p, i])) +(∫_298.15^T)(((Subscript[c, p])(T, p=1))/T)dT	  (2)

Data: I need  (Subscript[c, p])(T, p = 1atm) and s(T = 298.15 K, p = 1 atm). 

I know that (Subscript[c, p ])(T,p=1bar) for (Subscript[N, 2])is given by (Appendix 8.2):
 (Subscript[c, p])(T) = 6.449 + 1.413× (10^(-3)) T - 0.807 × (10^(-7))(T^2) (cal/Kmol)	(3)
In thermodynamic tables  I found:
	s(T = 298 K, p=1 atm) = 191.5 Joule/K mol		(4)
This is all I need to perform the calculation for the change in chemical potential

Mathematica plan: 
(1) Create a function cp[T] using Eq. (3)
(2) Create a function s[T] that returns 
	s(T, (Subscript[p, i]) ) =  s(T = 298.15 , (Subscript[p, i])) +(∫_298.15^T)(((Subscript[c, p])(T, p=1))/T)dT	
(3) Integrate  
	μ((Subscript[T, f]),(Subscript[p, i])) -μ((Subscript[T, i]),(Subscript[p, i])) = - (∫_(Subscript[T, i])^(Subscript[T, f]))s(T, (Subscript[p, i]))dT

Cell 1:  Get data for (Subscript[c, p])(T, 1atm) and s(298.15 K, 1 atm)

Define a function that give (Subscript[c, p]) (use Eq. 3 and the factor 4.184  to convert (Subscript[c, p]) from cal/mol K to J/mol K)
(Subscript[c, p])(T) = 6.449 + 1.413× (10^(-3)) T - 0.807 × (10^(-7))(T^2) (cal/K mol)

cp[T_]:=4.184*(6.449+1.413*10^-3*T-0.807*10^-7*T^2)(*Joule/K mol*)(*4.184 converts cal to Joule*)

cp[T] (*test that it works; gives (Subscript[c, p])(T, p=1atm) in J/K mol*)

Input the entropy at 298.15 K and p = 1 atm

s298=191.5 (*J/K mol*);

Cell 2: Define the function s[T] by using Eq. 2

s(T,1) = s(T=298, p=1) +(∫_298^T)(((Subscript[c, p])(T, p=1atm))/T)dT

Clear[a]
s[T_]=s298 +Integrate[cp[T]/T,{T,298,T},GenerateConditions->False] (*Joule/K mole*)

Check that the function works

s[x] (*check that it works; gives s(T,1atm) in J/mol K*)

Cell 3:  Calculate the change in μ, by using Eq. 1:
	Δμ  ≡  μ(T,p)-μ((Subscript[T, i]),p)  = - (∫_(Subscript[T, i])^(Subscript[T, f]))s(T, p)dT

I define a function 
	Δμ[T, p=1] = μ(T, p=1) - μ(298.15 K, p =1 atm) 
which gives me the change of the chemical potential as a function of temperature in the transformation
	{T=298.15, p = 1} -> {T, p = 1}.

Evaluate
	μ((Subscript[T, f]),(Subscript[p, i])) -μ((Subscript[T, i]),(Subscript[p, i])) = - (∫_(Subscript[T, i])^(Subscript[T, f]))s(T, (Subscript[p, i]))dT

Clear[Δμ]
Δμ[T_]=-Integrate[s[x],{x,298.15,T},Assumptions->T∈Reals&&T>0&&T<298]/1000. //Expand (*I divided by 1000 to get the change of μ in kJ/mol K*)

Mathematica explanation:  I used Assumptions to avoid getting a ConditionalExpression in the output.

Test that Δμ(298.15) is equal to zero

Δμ[298.15]

Plot[Δμ[T],{T,298,1000},AxesLabel->{"T","μ(T)-μ(298.15)"}]

The the chemical potential decreases with temperature meaning that μ(T, 1atm) is less than μ(298.15, 1 atm) if T < 298.15 K

Calculate μ(130 K, 1 atm) - μ(298.15 K, 1 atm):

Δμ[130]  (*kJ/mol*)

Cell 4: Calculate the same thing from μ = h - T s

I can use the equation 
	μ = h -T s
to calculate the change of μ:

μ(130K, 1atm) - μ(298 K, 1atm) = 
	(h(130 K, 1atm) - 130 s(130K, 1atm) - 
	(h(298 K, 1atm)- 298 s(298K,1atm))

Unfortunately F. Din, Thermodynamic Functions of Gases, vol. 3, Butterworth, London, 1961 gives only the values of h and s at 300K not at 298.15 K.  Therefore there should be a small difference between Δμ calculated above and the one obtained below
I have values for H at different temperatures
h(130K, 1atm) = 10.655 kJ/mol
h(300K, 1atm)=15.586 kJ/mol
s(130K, 1atm)=167.10/1000  kJ/mol K
s(300K, 1atm)=191.37/1000  kJ/mol K

Input the data:

h130=10.655;
h300=15.586;
s130=167.10/1000;
s300=191.37/1000; (*kJ for energy*)

Calculate μ(130 K, 1atm)-μ(300 K, 1 atm):

h130-130 s130-(h300 -300*s300)

This is very close to the value 30.4665 kJ/mol that I calculated by using 
	μ((Subscript[T, f]),p) -μ((Subscript[T, i]),p) = - (∫_(Subscript[T, i])^(Subscript[T, f]))s(T, p)dT	
for (Subscript[T, i]) =298.15 (while Din’s data is at 300 K)

WorkBookT.14.2. The change of chemical potential in an isothermal transformation.

Calculate the change in chemical potential of (Subscript[N, 2]) in the 
transformation 
{T = 130 K ,p=1 atm } ⟶{T = 130 K ,p=600 atm}

To evaluate the change of the chemical potential μ when the pressure changes from (Subscript[p, i]) to (Subscript[p, f]), while the temperature is held constant at T=(Subscript[T, f]), I use (see Eq. 14.12 in the textbook)
(Subscript[Δμ, 2])≡μ((Subscript[T, f]),(Subscript[p, f])) - μ((Subscript[T, f]),(Subscript[p, i])) = (∫_( v((Subscript[T, f]), (Subscript[p, i])))^( v((Subscript[T, f]),(Subscript[p, f]))))  v (Subscript[(((∂p)/(∂v))), T=(Subscript[T, f])])dv     (1)
Plan of computation:
1. I use the Benedict Webb Rubin (BWR) equation of state, since I am working at low temperature and high pressures.   The parameters in this equation, for (Subscript[N, 2]), are given in Appendix 6.  I chose the parameters that are good to a maximum pressure of 600 atm and low temperatures (second column in the table in Appendix 6, not the third column).
2. I'll use Mathematica to calculate the derivative in Eq. 1 and the integrand.
3. I'll calculate the molar volumes v((Subscript[T, i]) ,(Subscript[p, i])) and v((Subscript[T, f]),(Subscript[p, f])) by solving 		
	(Subscript[p, i]) = f((Subscript[T, i]), v)  to get (Subscript[v, i])    
and
	 (Subscript[p, f]) = f((Subscript[T, f]), v)  to get (Subscript[v, f])
 Here f is the BWR equation of state.
4. Then I ask Mathematica to perform the integral in Eq. (1)

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

Cell 1:  Get the necessary data.

The Benedict-Webb-Rubin equation of state (to calculate the pressure) is (this is taken form WB3HowToUseEqOfState

the Benedict-Webb-Rubin equation  (Eq. 2.11 in the textbook): 
	p=RTρ + (BRT-A-(ρ^2)(Subscript[C, 0])T-2  + (ρ^3)(bRT-a)+ a α (ρ^6)  + 
	(ρ^3)cT-2 (1 + γ (ρ^2)) exp[-((γρ)^2)]					(4)
 where 
 	 ρ =1/v.  
 Create a function that gives the pressure according to this equation

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

The parameters in the BWR equation for (Subscript[N, 2  ])are given in Appendix 6 (the second column has the right pressure and temperature range)

Clear[parametersBWR, A,B,C0,a, b, c, α,γ,R]
parametersBWR={A->1.1925,  B->0.0458,
			                                            C0->5.8891*10^3, a->0.0149,b->1.98154*10^-3,  c->5.48064*10^2,α->2.91545*10^-4,γ->7.5*10^-3, R->0.082 } ;

Check that the equation works (I use {} in the argument because the function expects a list; no harm is done if the list is empty):

pBWR[T,v,{}]

Put the parameter in and make sure that they are properly used:

pBWR[T,v,parametersBWR]

Calculate the pressure at room temperature when the molar volume is 24 liter/mol. The result should be close to 1.

pBWR[298,24,parametersBWR]

Cell 2:  Calculate the molar volumes in the initial state (130 K, 1 atm) and final state (130 K , 600 atm)

The function that does this is 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}] ];

Calculate the volumes at T=130 K and p=1 atm and at T=130 K and p=600 atm:

volumes={vi->volumeG[130,1,pBWR,parametersBWR],vf->volumeG[130,600,pBWR,parametersBWR]} (*liter/mol*)

I compare these with the results given in the p-v-t data by F. Din, Thermodynamic Functions of Gases, vol. 3, Butterworth, London, 1961:

	v(130 K, 1 atm) = 10.572 liter/mol
	v(130K, 600atm) = 0.036 liter/mol
	
The high pressure result is disappointing.  It indicates that the BWR equation of state does not work very well when T = 130 K and p = 600 atm.

Cell 3: Calculate the integral 
(Subscript[Δμ, 2])≡μ((Subscript[T, f]),(Subscript[p, f])) - μ((Subscript[T, f]),(Subscript[p, i])) =(∫_( v((Subscript[T, f]),(Subscript[p, i])))^( v((Subscript[T, f]),(Subscript[p, f]))))  v (Subscript[(((∂p)/(∂v))), T=(Subscript[T, f])])dv 
which is  Eq. (1)

ΔμWithPressure=(101.3325/1000) *Integrate[v D[pBWR[T,v,parametersBWR],v],{v,vi,vf}, GenerateConditions->False]

The factor 101.3325/1000 turns liter atmosphere  into kJ.  The analytic expression does not do us much good unless we happen to have the molar volumes in the initial and final state (from measurements). The numerical value of the integral is obtained by giving the volumes and the temperature:

ΔμWithPressure=ΔμWithPressure/.volumes/.T->130   (*kJ/mol*)

I can test whether this is in agreement with the data given in F. Din, Thermodynamic Functions of Gases, vol. 3,
Butterworth, London, 1961.  There I find:
	s(130 K, 600 atm) = 96.33 J/mol K. 
	h(130K,600 atm) = 7.305 kJ/mol
	s(130K, 1atm) = 167.10 J/mol K
	h(130K,1atm) = 10.655 kJ/mol

Input Din’s data all given at 130 K.

s1=167.10/1000 (*entropy at 1 atm and 130K in kJ/mol K*);
s600=96.33/1000(*entropy at 600 atm and 130K in kJ/mol K*);
h1=10.655 (*enthalpy at 1 atm and 130K in kJ/mol K*);
h600=7.305 (*enthalpy at 1 atm and 130K in kJ/mol K*);

Calculate the change 
	Δμ = μ(130 K, 600 atm, ) - μ(130 K, 600 atm)  
by using the formula 
	μ = h -Ts

Δμ=h600-130*s600-(h1 -130*s1)(*kJ/mol*)

The value of 5.85 kJ/mol that I got is  satisfactory.

WorkBookT.14.3. The change of chemical potential when both the temperature and the pressure are changing.

Calculate the change of chemical potential of (Subscript[N, 2]) when both pressure and temperature change.

the state changes from
	{(Subscript[T, i])=298.15 K, (Subscript[p, i])=1 atm}
to 
	{(Subscript[T, f])=130K and (Subscript[p, f]) = 600 atm}

In  WorkBook T.14.1 I calculated the change in chemical  potential when the temperature changes from 298.15 K to 130 K and the pressure was constant and equal to 1 atm.  The result was 

μ[130K, 1atm] - μ[298.15 K, 1 atm] = 30.4952 kJ/mol 			(1)

In WorkBookT.14.2 I calculated that

μ[130K, 600 atm] - μ[130 K, 1 atm] = 6.19453 kJ/mol			(2)

The first transformation is an isobaric transformation from initial to the final temperature.  The second is an isothermal transformation from the initial pressure to final pressure at constant final temperature.  

Adding Eq. 1 and Eq. 2 gives the answer we are seeking:

μ[90 K, 600 atm] - μ[298.15 K, 1 atm] = 30.47 + 8.01 = 38.48    kJ/mol

WorkBook T.14.4. Optional.  Another example of change of chemical potential in an isothermal transformation.

Calculate the change of chemical potential for (Subscript[CO, 2]) when the pressure changes from 1 atm to 100 atm and the temperature is held constant and equal to 298 K.  Use the Benedict Web Rubin equation. This calculation is similar to the one performed earlier in this file, but under new conditions for which we expect the equation of state to be more accurate.

The equation giving the change of chemical potential with pressure, in an isothermal transformation (T is held constant) is given by (Eq. 14.12 in textbook)
μ(T=298, p=100) - μ(T=298, p=1) = (∫_( (Subscript[v, i]))^( (Subscript[v, f]))) v (Subscript[(((∂p)/(∂v))), T])dv	(1)
The volumes (Subscript[v, i])and (Subscript[v, f])are calculated by solving the equations
	(Subscript[p, i])=f(298.15,(Subscript[v, i]))   							(2)
and 
	(Subscript[p, f])=f(298.15,(Subscript[v, f]))								(3)

Remove["Global`*"];

Cell 1:  The Benedict-Webb-Rubin equation is

p=RTρ + (BRT-A-(ρ^2)(Subscript[C, 0])T-2  + (ρ^3)(bRT-a)+ a α (ρ^6)  + 
	(ρ^3)cT-2 (1 + γ (ρ^2)) exp[-((γρ)^2)]	
 where  ρ =1/v.  
The constants appearing in these equations are given in Appendices  4-6 of the textbook

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

The parameters for the equation of state for (Subscript[CO, 2]) are:

Clear[param]
param ={A->2.51604,B->4.48842 10^(-2),C0->1.474405*10^5,a->1.3681*10^(-1),b->4.12381*10^(-3),c->1.49183*10^4,α->8.4658*10^(-5),γ->5.252*10^(-3),R->0.082}

Perform the usual test to find out whether the results are reasonable. For example, I happened to type incorrectly one parameters and the pressure for T=298 and v=24 was about 95 atm. Clearly something was wrong.

pBWR[T,v,{}]

pBWR[T,v,param]

pBWR[298,22,param]

Cell 2:  Calculate the change in chemical potential with Eq. 1 
μ(T=298, p=100) - μ(T=298, p=1) = (∫_(Subscript[v, i])^(Subscript[v, f]))v (Subscript[(((∂p)/(∂v))), T])dv

I need to know the molar volumes of the gas in the initial and the final state.  I calculate them by using the function volumeG provided in WB3

volumeG[T_,p_,f_,parameters_]:=
	Module[{ eq,sol},
Clear[C0,A,B,a,b,c,α,γ,R];
	eq=p==f[T,v ,parameters];

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

The initial volume is 24.44 liter/ mol and the final one is 0.133 liter/mol.  This is quite a large compression.

Calculate the change in chemical potential from
Δμ =μ(298K,100atm) - μ(298K,1atm) = 24.2147(∫_(Subscript[v, i])^(Subscript[v, f]))(Subscript[v(((∂p)/(∂v))), T])dv  cal/mol

The factor 24.2147 converts liter atm to calories.

int=v D[pBWR[298,v,param],v];
Δμ=24.2147*Integrate[int,{v,vi,vf}, GenerateConditions->False]/.volRules (*cal/mol*)

WorkBook T.14.5. The change of chemical potential for an isothermal transformation: use three equations of state: ideal gas, van der Waals and Beattie-Bridgeman

Calculate the numerical change in the chemical potential of a mole of oxygen in a transformation in which the state changes from 
	{ (Subscript[T, i])=298K,(Subscript[p, i])=1atm} 
to
	{(Subscript[T, f])=298K,(Subscript[p, f])=100 atm }

For the change of μ with pressure (at constant temperature) I use Eq. 14.13 of the textbook:
μ(T=298, p=100) - μ(T=298, p=1) = (∫_( (Subscript[v, i]))^( (Subscript[v, f])))v (Subscript[(((∂p)/(∂v))), T])dv		(1)
The volumes (Subscript[v, 1])and (Subscript[v, 2])are calculated by solving the equations
	(Subscript[p, i])=f(T=298,(Subscript[v, i]))   						(2)
and 
	(Subscript[p, f])=f(T=298,(Subscript[v, f]))							(3)
where p=f(T,v) is the equation of state of the system.

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

Cell 1:  The equation of state is the ideal gas equation.

Use Equation 14.21:

	μ((Subscript[T, i]), (Subscript[p, f])) -  μ((Subscript[T, i]), (Subscript[p, i])) =  R T  ln{(Subscript[p, f])/(Subscript[p, i])}			(4)
	
To get the result in calories use R =1.987  cal/mol K

ΔμId = 1.987 *298 Log[100./1] (*in cal/mol *)

The equation of state is the van der Waals equation.

The van der Waals equation is
	p= R T /(v-b) -a/(v^2)						(5)
To calculate the change in μ I use
	μ((Subscript[T, i]), (Subscript[p, f])) - μ((Subscript[T, i]), (Subscript[p, i])) =  n(∫_(Subscript[v, 1])^(Subscript[v, 2]))v (Subscript[(((∂p)/(∂v))), T])dv			(6)

Plan:  (1) Define a function that gives the pressure according to the van der Waals equation for oxygen. From Appendix 4.1 I get:
	a = 1.387 (liter^2)bar/(mol^2)   
	b = 0.03183 liter/mol
From Appendix 1,  1 bar = 0.987 atm. I need this to convert (liter^2)bar/(mol^2) to (liter^2)atm/(mol^2)

I calculate the volumes (Subscript[v, i])and (Subscript[v, f]) by solving the equations:
	1 = R 298.15 /((Subscript[v, i])-b) -a/((Subscript[v, i])^2)-b) 			(7)
and 
	100 = R 298.15 /((Subscript[v, f])-b) -a/((Subscript[v, f])^2)			(8)

Cell 2:  Get the van der Waals equation of state and the parameters for oxygen.

Input the parameters in the van der Waals equation, for oxygen

Clear[paramO2]
paramO2={a->1.387*0.987, b->0.03183,R->0.082}

Copy the function pvdw from WorkBookCh3

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

Check that you use the function pvdw correctly:

pvdw[T,v,{}]

pvdw[298,24,paramO2]

Cell 3:  Use the function volumeG (defined below) to calculate the initial and the final molar volumes (Subscript[v, i]) and (Subscript[v, f])

volumeG[T_,p_,f_,parameters_]:=
Module[{R=0.082 (*liter atm/mol*),sol},
	v/.FindRoot[ p==f[T,v ,parameters],{v,0.082 T/p}] ]

Test it

volumeG[298,1,pvdw,paramO2]

Define replacement rules for the volumes

volumes={vi->volumeG[298,1,pvdw,paramO2],
	vf->volumeG[298,100,pvdw,paramO2]}

Cell 4:  Calculate the change of chemical potential from Eq. 4.

The coefficient 24.2147 converts liter atmosphere into cal (see Appendix 7)

ΔμVdW=  24.2147* Integrate[ v D[pvdw[298,z,paramO2],z]/.z->v ,{v,vi,vf},GenerateConditions->False]/.volumes//Chop (*cal/mol*)

I used Chop to remove 0 I from the output.

Note the difference: if I use the ideal gas equation, the change in the chemical potential is 2726.84;  if I use the van der Waals equation it is 2672.26 cal/mol. This does not seem much.  But it makes a difference when we calculate equilibrium constants.

The equation of state is the Beattie - Bridgeman equation.

Cell 5:  Input the parameters for oxygen, in the Beattie-Bridgeman equation, from Appendix 5

paramBB={A->1.4911,a->0.02562,B->0.04624,
 b->0.004208, c->4.8*10^4}

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

Check that you use the equation correctly

pBB[298,22,paramBB]

Cell 6:   Calculate the initial and the final molar volumes:

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

vol={vi->volumeG[298,1,pBB,paramBB],
		vf->volumeG[298,100,pBB,paramBB]}

Cell 7:  Calculate the integral  (∫_(Subscript[v, i])^(Subscript[v, f]))v (Subscript[(((∂p)/(∂v))), T])dv;

The factor 24.2147 converts liter atm to calorie

ΔμBB=  24.2147* Integrate[ v D[pBB[298,z,paramBB],z]/.z->v ,{v,vi,vf},GenerateConditions->False]/.vol

This is pretty close to the van der Waals result which was 2672.26

WorkBookT.14.6. Optional.  Create a function that gives the temperature dependence of chemical potential when the pressure is constant

Remove["Global`*"];

Derive an equation giving the change of chemical potential of gaseous oxygen, when the temperature changes from 298 K to T and the pressure is held constant at p = 1 atm.  Plot   μ (T, 1 atm) - μ (298 K, 1 atm).

Plan: I will use the equation 14.6:
	μ((Subscript[T, f]),p) -μ((Subscript[T, i]),p)  = - (∫_(Subscript[T, i])^(Subscript[T, f]))s(T, p)dT			(1)
To calculate how s depends on temperature I use Eq. 14.9:
	s(T,p) =  s(T=298, p = 1) +(∫_298^T)(((Subscript[c, p])(T, p=1))/T)dT		(2)

Data: I need  (Subscript[c, p])(T, p=1atm) and s( T=298.15 K,p=1 atm). 
From the NIST data base I have  

(Subscript[c, p])(T) = 25.503 + 13.612 × (10^(-3)) T - 42.553 × (10^(-7))(T^2) (Joule/Kmol)		(3)
and
	s(T=298 K,p=1 atm) = 205.147 Joule/K mol			(4)

Mathematica plan: 
(1) Create a function cp[T] that returns Eq. (3)
(2) Create a function s[T] that returns s(T) calculated from Eq. 2.
(3) Integrate Eq. (1)

Cell 1:  Input data for (Subscript[c, p])(T,1atm) and s(298.15 K, 1 atm)

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

cp[T_]:=25.503 + 13.612  10^-3T - 42.553   10^-7  T^2  (*Joule/K mol*)

s298=205.147 (*J/K mol*);

Cell 2: Define the function s[T]:

s[T_]:=205.147 +Integrate[cp[y]/y,{y,298,T},GenerateConditions->False] (*Joule/K mole*)

s[T]

Cell 3:  Calculate the change in μ, from μ(T,p)-μ((Subscript[T, i]),p)  = - (∫_(Subscript[T, i])^(Subscript[T, f]))s(T, p)dT

I can define a function Δμ[T] which gives me the change of the chemical potential as a function of temperature

Clear[Δμ]
Δμ[T_]=-Integrate[s[x],{x,298,T}, GenerateConditions->False]/1000. //Expand (*kJ/mol K*)

Δμ[600]  (*kJ/mol*)

Cell 4: Plot the dependence of Δμ in T

Plot[Δμ[T] ,{T,298,2000}]

WorkBookT.14.7. Optional. Calculate the change of chemical potential of (Subscript[CH, 4]) in a transformation  from the initial state
 	   (Subscript[σ, i])={(Subscript[T, i])=298K, (Subscript[p, i])=1atm} 
to the final state
	 (Subscript[σ, f])={(Subscript[T, f])=400K, (Subscript[p, f])=100atm}.

Since the chemical potential is a function of state, its change in an arbitrary transformation is independent of path. This means that I can chose a convenient path for calculating the change.   I pick here the same path that I used when I calculated the change in entropy or enthalpy.  It consists of two succesive transformations:
	Transformation I is from
		 (Subscript[σ, i])={(Subscript[T, i])=298 K,(Subscript[p, i])=1 atm} 
to 
		(Subscript[σ, 1])= {(Subscript[p, i])=1 atm,(Subscript[T, i]) =400K} 
and it is at constant pressure (the path is p= 1 atm). 
	Transformation II is from
		 (Subscript[σ, 1])= {(Subscript[p, i])=1 atm,(Subscript[T, f]) =400 K} 
to 
		(Subscript[σ, f])={ (Subscript[p, f])=100 atm, (Subscript[T, f])=400K} 
and it is at constant temperature (the path is T=400 K)
	The total change of chemical potential is the sum of the change in transformation I and that in transformation  II.
	Implementation: The change in chemical potential, in the transformation consisting of the two segments described above is:
	Δμ ≡ μ((Subscript[T, f]),(Subscript[p, f]))-μ((Subscript[T, i]),(Subscript[p, i])) =  (Subscript[Δμ, I])+(Subscript[Δμ, II])			(1)
where 
	(Subscript[Δμ, I])= – (∫_( (Subscript[T, i]))^( (Subscript[T, f])))s(T,(Subscript[p, i]))dT					(2)
is the change of chemical potential in transformation I and 
	(Subscript[Δμ, II])= (∫_(v((Subscript[T, f]),(Subscript[p, i])))^( v((Subscript[T, f]),(Subscript[p, f]))))(Subscript[v(((∂p)/(∂v))), T])dv					(3)
is the change of chemical potential in the transformation II. Here v((Subscript[T, f]),(Subscript[p, i])) is the molar volume of the material at the temperature (Subscript[T, f])and the pressure (Subscript[p, i]). 
	To evaluate (Subscript[Δμ, II]) I need a good equation of state.  This can be used to calculate (Subscript[v(((∂p)/(∂v))), T]) and the molar volumes v((Subscript[T, f]),(Subscript[p, i])) and v((Subscript[T, f]),(Subscript[p, f])). Such molar volume calculations have been done in many of the previous chapters. For liquids or solids (Subscript[(((∂p)/(∂v))), T])is so small that (Subscript[Δμ, II]) is negligible (except when (Subscript[p, f]) -(Subscript[p, i]) is much larger than values normally used in chemistry).
	To evaluate (Subscript[Δμ, I]) I need to know how the entropy s(T,(Subscript[p, i])) depends on temperature at the pressure (Subscript[p, i]).  The NIST data base gives equations for the function s(T,1atm) for many compounds.  If you work with a compound for which such equation is not available you can calculate s(T,(Subscript[p, i])) from (see Chapter 11):
	s(T,(Subscript[p, i]))=s((Subscript[T, 0]),(Subscript[p, i])) +(∫_(Subscript[T, 0])^T)(((Subscript[c, p])(T,(Subscript[p, i])))/T)dT					(4)
In most practical cases (Subscript[p, i]) = 1 atm and (Subscript[T, 0]) = 298.15 K. Equations for (Subscript[c, p])(T,1atm) are given in the NIST data base and Appendix 8. The entropy at (Subscript[T, 0])=298 Kis  given in the NIST tables. 

Mathematica strategy: 
(1) Calculate (Subscript[Δμ, I])
(2)Calculate (Subscript[Δμ, II])
(3) Add them up to get Δμ
The strategy for implementing the first two steps is outline below

Remove["Global`*"];

Cell 1: Calculate the change of chemical potential (Subscript[Δμ, I]), in the transformation I.

Along path I the temperature changes from 298 K to 400 K, at constant pressure (Subscript[p, i])=1 atm. Therefore the entropy dependence on temperature is given by Eq. (4), which becomes:
s(T,p=1atm) = s(T=298.15K, p=1atm) + (∫_298^( 400))(((Subscript[c, p])(T, p=1atm))/T)dT	(5)

Plan of evaluation of Eq. (5):  
	(1) Write a function that gives the dependence of (Subscript[c, p]) on T at p=1 atm.  The heat capacity of (Subscript[CH, 4]) is (from Appendix 8.2):
	(Subscript[c, p])(T, p= 1 atm) = 4.171 + 14.450×(10^(-3))T + 1-0.267×(10^(-6))(T^2) cal/mol K       (6)
If T is in K the result for (Subscript[c, p]) is in cal/mol K
	(2) Calculate the integral in Eq. (5)
	(3) Add  s(T=298 K, p=1atm) = 188.66 J/mol K  to the integral (see Eq. 5)

Create the function cp[T] that returns the expression Eq. (6) and define the entropy at 298K

cp[T_] := 4.171 + (14.450/10^3)*T - (0.267*T^2)/10^6 (*cal/mol K*)
s298=188.66  (*Joule/mol K*);

Calculate s(T,p=1atm) from Eq. 5 and define a function s[T] that gives the entropy as a function of temperature.  The factor 4.184 converts the result from cal to Joule.

s[T_]:= s298+4.184*Integrate[cp[y]/y, {y,298,z},GenerateConditions->False]/.z->T (*J/mol K*)

Check that the function works properly

s[q]

Calculate the (Subscript[Δμ, I])from Eq. 2

ΔμI=-Integrate[s[T],{T,298.15,400}] (*Joule/K mol*)

Cell 2: Calculate the change of chemical potential in transformation II.

Plan for the calculation of (Subscript[Δμ, II]):
(1) To calculate the integral in 
	(Subscript[Δμ, II])= (∫_(v((Subscript[T, f]),(Subscript[p, i])))^( v((Subscript[T, f]),(Subscript[p, f]))))(Subscript[v(((∂p)/(∂v))), T])dv
which is Eq. (3),  I need the integration limits v((Subscript[T, f]),(Subscript[p, i])) and v((Subscript[T, f]),(Subscript[p, f])).  I calculate v((Subscript[T, f]),(Subscript[p, i])) by solving the equation (Subscript[p, i])=f((Subscript[T, f]), v) for v, where f is the Benedict Webb Rubin (BWR) equation of state (see Chapter 3 for examples of such calculations). v((Subscript[T, f]),(Subscript[p, f])) is obtained similarly, by solving (Subscript[p, f])=f((Subscript[T, f]), v) for v. The  function VolumeG  calculates the molar volume by using the prescription just described. 
(2)  I will use Mathematica to calculate the integrand v (Subscript[(((∂p)/(∂v))), T])from the BWR equation of state;
(3) Now I have all the information needed to  evaluate the integral Eq. (3).
	The data needed:  For these calculations I only need the BWR equation of state
	p =R T ρ + ((Subscript[B, 0])R T-(Subscript[A, 0])-(Subscript[C, 0])/(T^2))(ρ^2)+ (b R T -a) (ρ^3)+a α (ρ^6)+
	 ((cρ)^3)(1+γ(ρ^2))  exp[-((γρ)^2)]/(T^2)						(6)
Here ρ = 1/v.   For (Subscript[CH, 4])the parameters are:
(Subscript[A, 0])=1.855, (Subscript[B, 0])=0.0426, (Subscript[C, 0])=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, the equation (6) gives the pressure in atmospheres.

Copy the BWR equation from WorkBook3

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]

Input the parameters in the BWR equation for methane

param={A->1.79894, B->0.0454625,C0->3.18382*10^4,a->0.04352,
		b->0.00252033,c->3.5878*10^3,α->3.30*10^-4,γ->1.05*10^-2};

Check PBWR works well by using a high temperature and high molar volume, for which the ideal gas equation should work well (the package provides PIdeal[T,v] for calculating the pressure from the ideal gas law

pBWR[1000, 50,param]

The equation PBWR seems to work OK.

Use the function VolumeG, from the package, to calculate the molar volumes in the states {T=400 K, p=1 atm} and {T=400 K, p=100 atm}, which are the initial and the final states along path II.

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

Calculate the molar volumes

vol={vi->volumeG[ 400,1,pBWR,param],vf->volumeG[ 400,100,pBWR,param]}

Now calculate the integral 
	(Subscript[Δμ, II])= (∫_(v((Subscript[T, f]),(Subscript[p, i])))^( v((Subscript[T, f]),(Subscript[p, f]))))(Subscript[v(((∂p)/(∂v))), T])dv
(this is Eq. (3) ) to obtain (Subscript[Δμ, II])

ΔμII=4.184*101.325*Integrate[v D[pBWR[T, v,param], v]/. T -> 400, {v, vi, vf},GenerateConditions->False]/.vol (*converted to Joule*)

The factor 101.325 converts liter atm into caland 4.184 converts cal into Joule  (see Appendix 9) .

Cell 3: Calculate the total change in chemical potential Δμ= (Subscript[Δμ, I])+(Subscript[Δμ, II])

The total change in chemical potential, when the system is transformed from {T=298 K , p=1 atm} to {T=400K , p=100 atm}is the sum of the changes  on the two intermediate path

{ΔμI,ΔμII,Δμ=ΔμI+ΔμII}

The chemical potential changes quite a bit when we change pressure.  This has important consequences for applications to chemical equilibrium.

Cell 4:  As an exercise I calculate (Subscript[Δμ, II]) by using ideal gas law instead of the BWR

The function  v[T,p], defined below, gives the volume for an ideal gas

v[T_,p_]:=With[{R=0.08205},R  T/p]

volumes={vi->v[298.15,1],vf->v[400,100]}

The change in chemical potential in the transformation 
	(Subscript[σ, 1])={(Subscript[T, f])=1000K,p=1atm} ⟶ (Subscript[σ, f]){(Subscript[T, f])=1000K,(Subscript[p, f])=100atm}
is

ΔμIIIdeal=4.184*101.325*Integrate[v D[0.082  T/v,v] /. T -> 400, 
   {v, vi, vf},GenerateConditions->False]/.volumes (*Joule/K mol*)

The result obtained by using ideal gas law not close to that obtained by using a good equation of state:

{ΔμII,ΔμIIIdeal}

Had we worked at lower temperature or higher pressure, the error made by using ideal gas law would be worse

WorkBookT.14.8.  Find the equation for the fugacity of a gas that satisfies the van der Waals equation of state

Define the van der Waals equation

pvdw[T,v]:= R T/(v-b)- a/v^2

Use
	R T ln(f(T,p)/(Subscript[p, 0])) = (∫_(v(T,(Subscript[p, 0])))^( v(T,p)))v   (Subscript[(∂p/∂v ), T]) dv

v D[pvdw[T,v],v]

Calculate 
	R T ln(f(T,p)/(Subscript[p, 0])) = (∫_(v(T,(Subscript[p, 0])))^(v(T,p)))v   (Subscript[(∂p/∂v ), T]) dv

rtlnf = Integrate[v D[pvdw[T,v],v],{v,v0,v},GenerateConditions->False]

fTV= p0 Exp [rtlnf/(R T)]