© 2022 Horia Metiu

WorkBookT.6:  Heat

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

WorkBookT6.1  Fit the temperature dependence of (Subscript[c, p]) to a polynomial

Cell 1:  Input the data

At a pressure of one atm the heat capacity at constant pressure of (Subscript[CO, 2]) is given in the following table:

({{"t C", "\!\(C\_p\) Kcal/Kg"}, {-80, 0.183000000000000007`}, {-60, 0.18700000000000001`}, {-40, 0.191000000000000014`}, {-20, 0.195000000000000017`}, {0, 0.199000000000000021`}, {20, 0.203000000000000024`}, {40, 0.206999999999999984`}, {60, 0.210999999999999987`}, {80, 0.214999999999999991`}, {100, 0.218999999999999994`}, {120, 0.223}, {140, 0.226000000000000023`}})

• The data are from "Thermodynamic Functions of Gases" by F. Din.  Find a polynomial in T that fits this data well. 
• Weierstrass theorem (from calculus) says that any continuous function defined on a finite interval can be fitted, as well as desired, by a polynomial

Create a  list of temperatures in degree Kelvin

TData(*in degrees K*) =Table[i+273.16,{i,-80,150,20}]

Create the list cp containing the values of (Subscript[c, p ])corresponding to the temperatures in the list T

cpData(* in Kcal/kg K*)={0.183,0.187,0.191,0.195,0.199,0.203,
		0.207,0.211,0.215,0.219,0.223,0.226};

I want to use cal/mol K as units for (Subscript[c, p]).  I have a good reason for doing this:  I know the order of magnitude of (Subscript[c, p]) in cal/mole K which means that if I work in these units I can verify that my results are reasonable. 
 
To convert the cpData from Kcal/Kg K to cal/mole K I have to multiply (Subscript[c, p]) with the molecular mass of (Subscript[CO, 2])

cpData = cpData *(12+2 *16)                              (*cal/mol K*)

These values have the right order of magnitude.

I want to use Mathematica function Fit  to generate a formula that give (Subscript[c, p]) as a function of T. The way the data is stored (see the line below) is not appropriate for use in Fit.

{TData,cpData}

Create a data list that has the structure needed for use in the Mathematica function Fit

TcpData=Transpose[{TData,cpData}]

Cell 2:  Use Mathematica function Fit to fit the data with a polynomial

I plan to use Mathematica function Fit to fit, with a polynomial, the dependence of (Subscript[c, p]) on temperature. The syntax is 
	Fit[{ {(Subscript[T, 1]),(Subscript[c, p])((Subscript[T, 1]))},...,{(Subscript[T, n,])(Subscript[c, p])((Subscript[T, n]))}} , {1, x, x^2,...,x^m}, x]
The data is the list in the curly blue bracket.  This consists of lists with two elements the temperature and the corresponding heat capacity.   The list rendered in purple tells the function Fit to fit the data with 
	f(x) = (Subscript[a, 0]) + (Subscript[a, 1]) x +  ... + (Subscript[a, m]) (x^m) 
Fit will determine the magnitude of the coefficients (Subscript[a, i]) so that the function f(T) fits (Subscript[c, p])(T )
	(Subscript[c, p])(T ) ≃ f(T)

Define the function monomials[x][m] that creates the list of monomials {1, x, x^2,...,x^m} used by the function Fit.

Clear[monomials]
monomials[x_][m_]:=Module[{n},tb=Table[x^n,{n,0,m}]]

Test that it works: generate all monomials up to and including (x^6)

monomials[x][6]

Use Mathematica function Fit to define function cpFit[m] that approximates the data given in TCpData with a polynomial of order m. The syntax is 
Fit[{ {(Subscript[T, 1]),(Subscript[c, p])((Subscript[T, 1]))},...,{(Subscript[T, n,])(Subscript[c, p])((Subscript[T, n]))}} , {1, T, (T^2),(T^3),...,(T^m)}, x]

Clear[cpFit]
cpFit[x_ ][m_Integer]:= Fit[ TcpData,monomials[T][m],T]

cpFit[T][6]

Look at a fit of the data in the list  CpTData with a six order polynomial

cpCalculated=cpFit[ T][6] /.T->#&/@TData

cpFit[300][10]

Cell 3:  Test whether the fit is good.

There are two way of testing whether the fit is good: 
• make a graph with the data and the fit
• make a table with the data, the fit and the error

Make a table with the experimental values of (Subscript[c, p]), the values calculated with the equation I derived, and the error made by using my equation:

TableForm[MapThread[{#1,#2,#3,(#2-#3)*100/#3}&,{TData,cpData,cpCalculated}],
TableHeadings->{None,{"T K","(Subscript[c, p]) exp","(Subscript[c, p])calc", "% error"}}]

The errors are small.  The equation we have obtained is fairly good.

Make a plot with the data by using ListPlot

listPl=ListPlot[TcpData,PlotStyle->{RGBColor[1,0,0],PointSize[0.025]},AxesLabel->{"T K","(Subscript[c, p]) cal/mole K"}] ;

Plot the fitting function

plfit=Plot[cpFit[T][10]//Evaluate,{T,193,420},PlotStyle->RGBColor[0,0,0]];

Show[{listPl,plfit}]

Cell 4:  Compare my fit with a fit found in the literature

I found in the literature the equation 
	(Subscript[c, p]) = 5.152 + 15.224×(10^(-3))T – 96.81× (10^(-7))(T^2) + 2.313×(10^(-9))(T^3)
Here T is in degrees K and the result for (Subscript[c, p ])is in cal/mol K. Let us compare this equation to ours. Define a function based on the equation from the literature:

cpLit[T_] := 5.152 + (15.224*T)/(10^3) - (96.81 T^2)/10^7 + (2.313/10^9)*T^3

cpLit[x]

This function gives the following values for the specific heat:

cpFromLit=cpLit/@TData

Compare this with the experimental values stored in cp:

TableForm[MapThread[{#1,#2,#3,(#2-#3)*100/#3}&,{TData,cpData,cpFromLit}],
TableHeadings->{None,{"T K","(Subscript[c, p]) exp","(Subscript[c, p])FromLit", "% error"}}]

Overall, my function does a better job then the one I found in the literature. This is not surprising. The equation from the literature fitted data in the range 300 K to 1500 K. I fitted  data taken in the temperature range for which I am testing the equation.

WorkBookT6.2  Fit the pressure dependence of (Subscript[c, p]) to a polynomial

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

I want to fit the pressure dependence of (Subscript[c, p])with the polynomial
	(Subscript[c, p])(p) = a + b p + c (p^2) + d (p^3) + etc.
The pressure dependence of the heat capacity is important for ammonia synthesis and for methanol synthesis because we need to calculate the heat of reaction and the equilibrium concentrations at high pressure.

Cell 1:  Data on (Subscript[c, p]) versus pressure when T = 330 K

I pick data from Din’s book and make a list of the form
	{ {(Subscript[p, 1]), (Subscript[c, p])((Subscript[p, 1]))},{(Subscript[p, 2]),(Subscript[c, p])((Subscript[p, 2]))},...}

pCpData={{1,13.62},{2,13.71},{3,13.80},{4,13.89},{5,13.97},{6,14.01},{7,14.12},{8,14.21},{9,14.30},{10,14.40}{15,14.93},{20,15.57},{25,16.31},{30,17.17},{35,18.25},{40,19.73},{45,21.85},{50,24.92},{55,29.96},{60,38.40},{65,46.45},{70,48.54},{75,46.80},{80,43.64},{85,40.35},{90,37.40},{95,34.94},{100,32.92}};

Let’s make a plot the data

ListPlot[pCpData,AxesLabel->{"p atm","(Subscript[c, p]) cal/mole K"},PlotLabel->"The dependence of (Subscript[c, p]) on pressure\n"]

I want to make a list with the pressures  contained in  the list pCpData.  I use pattern matching and I call the list of pressures pListData

pData=pCpData/.{x_,y_}->x

I want to make a list with the heat capacities

cpData=pCpData/.{x_,y_}->y

Cell 2:  Fit the data with a polynomial

I plan to use Mathematica function Fit to fit, with a polynomial, the dependence of (Subscript[c, p]) on pressure. The syntax is 
	Fit[{ {(Subscript[p, 1]),(Subscript[c, p])((Subscript[p, 1]))},...,{(Subscript[p, n]),(Subscript[c, p])((Subscript[p, n]))}} ,{1, x, x^2,...,(x^m)}, x]
The data is in red and it is a list of lists.  The list rendered in purple tells the function Fit to fit the data with 
	f(x) = (Subscript[a, 0]) + (Subscript[a, 1]) x +  ... + (Subscript[a, m]) (x^m) 
Fit will determine the coefficients (Subscript[a, i]) so that the function f(p) fits (Subscript[c, p])(p)
	(Subscript[c, p])(p) ≃ f(p)

Define the function monomials[n] that creates the monomials used in the function Fit.

Clear[monomials]
monomials[x_][n_]:= Table[x^m,{m,0,n}]

Test that it works: generate all monomials up to (x^6)

monomials[p][6]

Use Mathematica function Fit to define function cpFit[n] that approximates the data given in pCpData with a polynomial of order n.

Fit[pCpData,monomials[p][2],p]

Clear[cpFit]
cpFit[p_][n_Integer]:= Fit[pCpData,monomials[x][n],x] /.x->p

Look at a fit of the data in the list pCpData with a six order polynomial

cpFit[p][6]

This is a six order polynomial fit.  If I make p = 10,  I get the heat capacity  (Subscript[c, p]) given by the fit when p = 10 atm

cpFit[10][6]

If I look at the data in the list pCpData I find  {10, 14.40} which tells me that at a pressure of 10 atm,  (Subscript[c, p]) = 14.40 cal/mole K.  The fit with a six order polynomial gives (Subscript[c, p]) = 15.4539 cal/mole K.  I will explore below fits with polynomials of different orders.

I calculate (Subscript[c, p]) at pressure of 10 atm, fitted with polynomials of order 4, 6, 8, 10, 12, 14, 16.

Table[{n,cpFit[10][n] },{n,4,16,2}]

Calculate the error made by the fit:
 	 percentError = (((Subscript[c, p])(p) -(Subscript[c, p,fit])(p))/((Subscript[c, p])(p))) 100

I had to use a polynomial of order 10 to get a reasonable fit.

The fit with a third order polynomial is bad

percentError=100*(cpData -(cpFit[#][3]&/@pData))/cpData

The fit with a polynomial of order ten is fine.

percentError=100*(cpData -(cpFit[#][10]&/@pData))/cpData

WorkBookT6.3  The heat exchanged in an isobaric transformation.

Cell 1. Formulate the problem and collect the data you need.

How much heat is needed for changing the temperature of 100 (m^3) of (Subscript[CO, 2]) from 400 K to 1000 K,  if the pressure is held constant and equal to 3 atm ?   
	(Subscript[q, p]) = n(∫_( (Subscript[T, 1]))^( (Subscript[T, 2]))) (Subscript[c, p])(p=3 atm, T)dT								(1)
While heating the gas, the pressure is held constant, therefore I use (Subscript[c, p]) in this equation

In Appendix 8.2 I find that at a pressure of one atmosphere the (Subscript[c, p]) for (Subscript[CO, 2]) is given by
	 (Subscript[c, p]) = 5.152 +15.224 ×(10^(-3))T - 96.81 × (10^(-7))(T^2) + 2.313×(10^(-9))(T^3) 	(2)
 I assume that this equation is correct at a pressure of 3 atm.  I know that this is not quite right but I don’t have data at p = 3 atm.

To evaluate (Subscript[q, p]) with Eq. 1 I need to calculate the number of moles n.  n  does not change during the transformation (no chemical reaction, no leaks from the container).   This means that I can calculate n from the data on  the initial state, which is  
	V = 100 (m^3) , p = 3 atm and T = 400 K.   
To find n I use the equation of state p = f(T, v) to calculate the molar volume by solving the equation
	3 atm = f(400, v) 
Once I know v I calculate the number of moles n from 
	n = (V/v)  
To perform this calculation I need to pick an equation of state.  I will use the van der Waals equation 
	p = ((R T)/(v-b))- (a/(v^2))
For (Subscript[CO, 2]) 
	a = 3.592 (liter^2) atm/(mole^2)
and 
	b = 0.04267 liter/mole

Cell 2. Define (Subscript[c, p])(T) and the pressure according to van der Waals equation

cp=.
cp[T_]:= 
   5.152 +  15.224 10^(-3) T -  96.81 10^(-7) T^2  +  2.313  10^(-9) T^3

Check if this works properly:

cp[x]

Define a function that gives 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];

Define the parameters a and b for (Subscript[CO, 2]) (use Appendix 4)

vdwParameters ={a->3.592,b->0.04267}

Test that these equations work  properly:

Clear[t,v]
{pvdw[T,v,{}],pvdw[T,v,{a->3.592,b->0.04267}]}

Cell 3. Calculate the number of moles

Calculate the molar volume of the initial state by solving
	3 atm = pvdw[400, v]
You know that  V = 100 (m^3) in the initial state and you calculate the number of moles from
	n = (V/v)

To calculate the molar volume use the function volumeG which we defined in WBT3 and which I give 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}] ];

Calculate the molar volume in the initial state: p = 1 atm and T = 400 K

Clear[v]
rulev=v->volumeG[400,3,pvdw,{a->3.493,b->0.042}] (*liter/mol*)

The number of moles is V/v. Since v is given in liter/mole, I must use V in liter. I know that 
	V = 100 (m^3)=100×(10^3) liters. 
The number of moles is

nbMoles=100*10^3/v  /.rulev

Cell 4. Calculate the heat needed for changing the temperature of the gas

Use 
	(Subscript[q, p]) = n (∫_( 400)^1000)(Subscript[c, p])(p=3, T)dT

qp=nbMoles Integrate[cp[T],{T,400,1000}] (*calories*)

Summary:  One needs 5.46 (10^7) calories to heat 100 (m^3) of (Subscript[CO, 2]) from 400 K to  1000 K, while keeping the pressure constant and equal to 3 atm,

WorkBookT6.4. The heat produced when one mole of (Subscript[CO, 2]) is compressed from p = 1 atm to p = 100 atm at constant temperature  T = 350 K.

This is the kind of calculation that one needs to perform when designing a compressor. If I want to keep the temperature constant during compression I need to know how much heat I need to remove.

•  When I explain the applications of the second law I will show that the heat (Subscript[q, T]) generated by isothermal compression is given by
	(Subscript[q, T]) = T (∫_( (Subscript[v, i]))^( (Subscript[v, f]))) (Subscript[(((∂p)/(∂T))), v])dv										(1)
• This is a very handy equation because the derivative of p with T (at constant v) is easy to calculate from the equation of state. This is a lot simpler than doing calorimetric measurements. 
• A minor nuisance is that Eq. 1 will give the heat in liter atm. I want to get the heat in Joules. I have (see Appendix 7)  
	1 liter atm = 101.325 Joule
and
	(Subscript[q, T]) = 101.325 T (∫_( (Subscript[v, i]))^( (Subscript[v, f]))) (Subscript[(((∂p)/(∂T))), v])dv								(2)
Here p is in atmospheres, v in liter/mole and (Subscript[q, T]) is in Joule

Define the van der Waals equation of state for (Subscript[CO, 2])

pvdw[T_,v_,parameters_]:=With[{R=0.082 (*liter atm/mole K*)},
R T/(v-b)-a/v^2/.parameters]

The parameters for (Subscript[CO, 2]) are (Appendix 4)

vdwParameters={a->3.593,b->0.042}

Test that these equations work properly

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

Calculate the heat in Joule.  Use 
	(Subscript[q, T]) = 101.325 T (∫_( (Subscript[v, i]))^( (Subscript[v, f]))) (Subscript[(((∂p)/(∂T))), v])dv           Joule
Use the van der Waals equation to calculate  
	(Subscript[(∂p/∂T), v])

To make the code clear I calculate the derivative first and then I calculate the integral

derivative[T_,v_]:=D[pvdw[T,v,vdwParameters],T]

Clear[v,vi,vf]
heat=101.325  T Integrate[ derivative[T,v],{v,vi,vf}, GenerateConditions->False]

Use the function VolumeG to calculate the molar volume of the gas:

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

Clear[vf,vi]
ruleMolarVolumes={vi->volumeG[350, 1 (*atm*), pvdw,{a->3.593,b->0.042}],vf->volumeG[350, 100(*atm*), pvdw,{a->3.593,b->0.042}]}

The heat produced in an isothermal compression from 1 atm to 100 atm, at a constant temperature of T = 350 K, is

heat=heat/.ruleMolarVolumes/.T->350

The calculation using the ideal gas gave -13,400 Joule/mole. The percent error made by using ideal gas law versus the van der Waals equation is

(heat-(-13400))*100/heat

The error made by using the ideal gas law is 14.84 % .  If you underestimate your budget for coolant (you are likely to want to remove the heat produced by compression) by 14%, you are in trouble.  You can not use the ideal gas law for this kind of calculation (the final pressure is too high).  The van der Waals equation is not safe either at 100 atm. You should use the BWR equation.

WorkBookT.6.5.Calculate the temperature reached when heat is transferred from  a hot gas to a cold one, at constant volume.

I put in contact two containers, one having inside 
	1 (m^3)  of (Subscript[CO, 2]) at 1000 K and p = 1 atm, 
and the other having 
	2 m3of N2 at 298 K and p = 1 atm.  
Heat is transferred between them and no heat is allowed to escape into the environment.  The total volume is constant.  What temperature do the gases have when they reach thermal equilibrium? How much heat do they exchange?

• As I explained in Chapter 6.3, the temperature reached by the gases is calculated from the condition that the heat lost by the cooling gas is equal to the heat taken by the gas being warmed up. 
• I emphasize that this is true only if the system does no work. This is the case here because the volumes of the two containers are constant.  
• Therefore I can write 
	- n((Subscript[CO, 2]))(∫_1000^T)(Subscript[c, v])((Subscript[CO, 2]))dT = n((Subscript[N, 2]))(∫_298^T)(Subscript[c, v])((Subscript[N, 2]))dT 				(1)

and solve this equation to get the equilibrium temperature T. 
•  I must use (Subscript[c, v]), not (Subscript[c, p]), because the process takes place at constant volume.  
• To keep life simple I use the ideal gas law 
	p v = R T 
• I assume that the heat capacity is not pressure dependent . You can figure out, as an exercise, how the calculation will change if I were to use Benedict-Webb-Rubin equation.

Often I find data for (Subscript[c, p]) and I calculate (Subscript[c, v]) from
	(Subscript[c, v])=(Subscript[c, p]) -R										 		 (2)
which is valid for ideal gases. Since p = 1 atm and the temperature is fairly high, this is not a bad approximation.

In Appendix 8.2  I found:
	(Subscript[c, p])((Subscript[N, 2])) = 6.449 + 1.413×(10^(-3)) T - 0.0807×(10^(-6))(T^2)	cal/mol K	 	(3)
	(Subscript[c, p])((Subscript[CO, 2])) = 5.152 + 15.224×(10^(-3)) T - 9.681×(10^(-6))(T^2)	 cal/mol K	(4)
To obtain (Subscript[c, v]) I use Eq. 2 together with Eqs. 3 and 4:	
	(Subscript[c, v])((Subscript[CO, 2])) = (Subscript[c, p])((Subscript[CO, 2])) - 1.987   		cal/mol K      				(5)
	(Subscript[c, v])((Subscript[N, 2])) =   (Subscript[c, p])((Subscript[N, 2]))  -   1.987   		cal /mol K				(6)
In the equations above, I take R = 1.987   cal/K mol because these are the units used for (Subscript[c, p]).

To calculate the number of moles n (which are needed in Eq. 1), I need to find the molar volume v and use n = V/v, where V is the total volume of the gas.   I assume that the equation of state for both gases is the ideal gas law. Obviously, I am sacrificing accuracy to reduce the labor required for solving the problem. My excuse is that I want illustrate the principle of the calculations without being distracted by algebraic difficulties. The trouble is that in real life, if you ignore these complications, you get fired.

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

Cell 1: Create a function that generates Eqs. 5 and 6 for (Subscript[c, v]).

Create functions cvCO2[T] and cvN2[T] that, given the temperature, calculate the heat capacities according to Eqs. 3 and 4

Define the function cvN2[T] which calculates
	(Subscript[c, v])((Subscript[N, 2])) = 6.449 + 1.413×(10^(-3))T - 0.0807×(10^(-6))(T^2)- 1.987	cal/mol K	 	(3)

Clear[cvN2];
cvN2[T_]:=With[{R=1.987 (*cal/mol K*)}, 6.449+1.413*10^-3 T-0.0807*10^-6  T^2 - R ];
cvN2[x]

Define the function cvCO2[T] which calculates
	(Subscript[c, p])((Subscript[CO, 2])) = 5.152 + 15.224×(10^(-3))T - 9.681×(10^(-6))(T^2)- 1.987	cal/mol K		(4)

cvCO2[T_]:=With[{R=1.987 (*cal/mol K*)}, 5.152+15.224×10^-3 T-9.681×10^-6 T^2- R ];
cvCO2[x]

Cell 2:  Calculate the number of moles of gas in each reservoir.

It is in this step that using the ideal gas equation of state makes life easy.  The number of moles are given by 
	 n =  p V/ R T.

For (Subscript[CO, 2 ])I know that 
	V = 1(m^3) = 1000 liter, p = 1 atm, T=1000 K .  
This gives for the number of moles (use R = 0.082 liter atm/mole K)

nCO2=1(*p in atm*) 10^3(*V in liter*)/(0.082*1000 (*T inK*))

For (Subscript[N, 2]) I have
	 V = 2 (m^3) = 2000 liter, T = 298 K.  
The number of moles is

nN2=2 *10^3/(0.082*298)

Cell 3:  Calculate the equilibrium temperature

I must solve for T the equation 
	-(Subscript[n, (Subscript[CO, 2])])×(∫_1000^T)(Subscript[c, V])((Subscript[CO, 2]))dT= (Subscript[n, (Subscript[N, 2])])×(∫_400^T)(Subscript[c, V])((Subscript[N, 2]))dT		(1)

Calculate first the heat given away by (Subscript[CO, 2 ])when it goes from 1000 K to a temperature T:

Clear[T];
heatGiven=nCO2*Integrate[cvCO2[t],{t,1000,T}]

Calculate the heat received by (Subscript[N, 2]) when it goes from 298 K to T

heatReceived= nN2*Integrate[cvN2[t],{t,298,T}]

I can find the temperature reached by the system from the condition that the amount of heat received by N2 must be equal to the amount given away by (Subscript[CO, 2]):

eq= (heatGiven//Expand)+(heatReceived//Expand)==0

Remember that heat received from environment has positive  sign and heat given (to the environment) has a positive sign.    Solve the equation

sol=NSolve[eq,T]//Flatten

The only meaningful temperature is T = 445  K,  because the final temperature must be between the initial temperature (298 K) of the cold gas and the initial temperature of the hot gas (1000K).