Notebook 3.2: Predicting Boiling Points from Liquid/Vapor Gibbs Free Energy Functions

January 26, 2010

This notebook shows how Shomate functions for a liquid and vapor can be used to predict boiling points. It also introduces function programming and for loops.

* Download and install the Euler Math Toolbox.
* Cut and paste the code below into a plain text file.
* Save the file with an .en file extension.
* Double-click the file. It should open up in the Toolbox.

All of the notebooks in this series are specifically keyed to Atkins’ Physical Chemistry, 8th edition. Italics mark the items students had to fill in themselves.

----------------------------------------------------snip snip----------------------------------------------------------
 Notebook 3.2: Variation of the Gibbs Free Energy with Temperature
% In this notebook, we'll predict the boiling points of a liquid by looking for a 
% temperature that makes the Gibbs free energies of the vapor and liquid equal. 
% To compute the Gibbs free energy at standard state (constant pressure), 
% we'll integrate the equation
%               dG = -SdT
% We'll need to find S as a function of T to do the integral. A very accurate 
% empirical fit called the Shomate equation can be used here:
% S = A*log(t) + B*t + C*t^2/2 + D*t^3/3 - E/(2*t^2) + G;
% where we can find look up the Shomate coefficients A, B, C, D, E, and G on 
% the NIST Chemistry Webbook (http://webbook.nist.gov/chemistry). 
% When we're doing a complicated calculation that will be repeated many times, 
% it's convenient to define the calculation as a "function". In Euler, we can define
% a function with a series of statements like
%      function f(x)
%          ....
%      return y
%      endfunction
% where f is the function name, x is value passed to the function, and y is the
% calculated value of f(x) to be returned. For example, here is a function that 
% computes the pressure of an ideal gas from its temperature and molar volume:
>function Pideal(V,T)
$return R*T/V;
% We can compute the pressure of an ideal gas with V = 22.4 L/mol and 
% T = 273.15 K in atmospheres as
>Pideal(22.4, 273.15)
% which is exactly the same as typing
% Let's define a function that computes the entropy of steam at temperature 
% T using the Shomate equation. The Shomate coefficients for water and steam
% are taken from 
% http://webbook.nist.gov/cgi/cbook.cgi?ID=C7732185&Units=SI&Mask=7.
>function EntropySteam(T) 
$return A*log(t) + B*t + C*t^2/2 + D*t^3/3 - E/(2*t^2) + G;
% This function computes the entropy in SI units (J mol^-1 K^-1 ) 
% To check to see that we've typed all the coefficients in correctly, we can 
% compare our entropy at 298.15 with the CODATA Review value, 
% 188.835 +/- 0.010 J mol^-1 K^-1 
% (Cox, J. D.; Wagman, D. D.; Medvedev, V. A., CODATA Key Values for 
% Thermodynamics, Hemisphere Publishing Corp., New York, 1984, 1). 
% Looks good. Now let's calculate the entropy of liquid water, again using 
% Shomate equation coefficients from the NIST Webbook:
>function EntropyWater(T) 
$t = T/1000;
$A = -203.6060;
$B = 1523.290;
$C = -3196.413;
$D = 2474.455;
$E = 3.855326;
$G = -488.7163;
$return A*log(t) + B*t + C*t^2/2 + D*t^3/3 - E/(2*t^2) + G;
% Again, we see that the function correctly predicts the CODATA Review value
%  of the standard entropy of water at
% 298.15 K, which should be 69.95 ± 0.03 J mol^-1 K^-1.
% To integrate -S dT from 298.15 K to 500 K, we can type
% Typing "EntropySteam(T)" won't work; because Euler doesn't know we want to 
% integrate with respect to T. 
% We have to call the integration variable "x".  
% The integral gives us the change in Gibbs free energy between 298.15 and 
% 500 K. To get the Gibbs free energy itself at some temperature T, we can 
% add the Gibbs free energy of formation at 298.15 to the integral of -S dT
% from 298.15 K to T:
>function GSteam(T)
$DeltaGfo = -228570;
$return DeltaGfo - integrate("EntropySteam(x)",298.15,T)
% At 500 K, the Gibbs free energy of steam is
% Now let's build an array of Gibbs free energies for the gas (call it Ggas)
% and a corresponding array of temperatures. Suppose we want 200 points:
>Ggas = 1:200; T = 1:200;
% Unfortunately, we can't just type Ggas = GSteam(T), because 
% our GSteam function only calculates one Gibbs free energy at a time. We
% can use a "for" loop to process each individual element of the T and Ggas 
% arrays. The syntax is
%     for index=start to finish;  ... end
% For example, if we wanted to set up 200 temperature points, starting from
% 299.15 K and ending with 499.15 K, we could type
>for i=1 to 200; T[i] = 298.15 + i; end
% Now set up the Ggas array. We'd like the Gibbs free energy to
% be displayed in kilojoules per mole.
>for i=1 to 200; Ggas[i]= GSteam(T[i])/1000; end
% Plot the Gibbs free energy of the gas as a function of T:
>xlabel("Temperature (K)");
>ylabel("Gibbs Free Energy (kJ/mol)");
% Now let's add a liquid curve to the plot. First, define a function
% to calculate the Gibbs free energy of water at a specific temperature:
>function GWater(T)
$DeltaGo = -237130;
$return DeltaGo -integrate("EntropyWater(x)",298.15,T)
% Now set up an array to hold the Gibbs free energy data for the liquid
% (call it Gliq). 
>Gliq = 1:200;
>for i=1 to 200; Gliq[i] = GWater(T[i])/1000; end
% Add the liquid curve to the graph.
>title("Determination of Boiling Point from Vapor/Liquid Gibbs Free Energies");
% When the Gibbs free energies of the liquid and gas are equal,
% the liquid boils. To get the boiling point Tbp, then, we need
% to solve the equation GSteam(Tbp) = GWater(Tbp) for Tbp.
% In a previous notebook, we used Euler's root() function to solve
% equations like this. We need an initial guess for Tbp. Anything
% that isn't too far off will do, so we'll use the midpoint of our
% temperature range as a rough guess:
>bp = mean(T)
% Now use root to find the intersection of the two curves. 
>bp = root("GSteam(bp)-GWater(bp)",bp)
% ...which is only about .05 K different from the true standard
% boiling point of water.
% We can double-check by printing the Gibbs free energies of the
% liquid and gas at this temperature. They must be equal.
% Finally, let's mark the boiling point on the plot. We want to
% drop a line from the intersection at (bp, GSteam(bp)/1000) to
% the x axis at (bp, min(Ggas)). 
>label("Boiling Point",bp,Ggas[190]);

% Make a NEW notebook (using this one as a model) that locates the
% boiling point of some other compound by plotting the Gibbs free
% energy for the phases in equilibrium. Your notebook should include
% a graph set up and labeled like the one above. It should also
% include a numerical estimate of the boiling or melting temperature.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Google photo

You are commenting using your Google account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s

%d bloggers like this: