Gibbs energy minimization and the NIST webbook

December 26, 2011 at 01:10 AM | categories: optimization | View Comments

Gibbs energy minimization and the NIST webbook

Gibbs energy minimization and the NIST webbook

John Kitchin

Contents

In Post 1536 we used the NIST webbook to compute a temperature dependent Gibbs energy of reaction, and then used a reaction extent variable to compute the equilibrium concentrations of each species for the water gas shift reaction.

Today, we look at the direct minimization of the Gibbs free energy of the species, with no assumptions about stoichiometry of reactions. We only apply the constraint of conservation of atoms. We use the NIST Webbook to provide the data for the Gibbs energy of each species.

As a reminder we consider equilibrium between the species $CO$, $H_2O$, $CO_2$ and $H_2$, at 1000K, and 10 atm total pressure with an initial equimolar molar flow rate of $CO$ and $H_2O$.

Links to data

H2

H2O

CO

CO2

function nist_wb_constrained_min

T = 1000; %K
R = 8.314e-3; % kJ/mol/K

P = 10; % atm, this is the total pressure in the reactor
Po = 1; % atm, this is the standard state pressure

We define species like this.

We are going to store all the data and calculations in vectors, so we need to assign each position in the vector to a species. Here are the definitions we use in this work.

1  CO
2  H2O
3  CO2
4  H2
species = {'CO' 'H2O' 'CO2' 'H2'};

Heats of formation at 298.15 K

Hf298 = [
    -110.53  % CO
    -241.826 % H2O
    -393.51  % CO2
       0.0]; % H2

Shomate parameters for each species

   A             B           C          D          E       F          G       H
WB = [25.56759  6.096130     4.054656  -2.671301  0.131021 -118.0089 227.3665   -110.5271  % CO
      30.09200  6.832514     6.793435  -2.534480  0.082139 -250.8810 223.3967   -241.8264    % H2O
      24.99735  55.18696   -33.69137    7.948387 -0.136638 -403.6075 228.2431   -393.5224   % CO2
      33.066178 -11.363417  11.432816  -2.772874 -0.158558 -9.980797 172.707974    0.0]; % H2

Shomate equations

t = T/1000;
T_H = [t; t^2/2; t^3/3; t^4/4; -1/t; 1; 0; -1];
T_S = [log(t); t; t^2/2; t^3/3; -1/(2*t^2); 0; 1; 0];

H = WB*T_H; % (H - H_298.15) kJ/mol
S = WB*T_S/1000; % absolute entropy kJ/mol/K

Gjo = Hf298 + H - T*S; % Gibbs energy of each component at 1000 K

Gibbs energy of mixture

this function assumes the ideal gas law for the activities

    function G = func(nj)
        Enj = sum(nj);
        G = sum(nj.*(Gjo'/R/T + log(nj/Enj*P/Po)));
    end

Elemental conservation constraints

We impose the constraint that all atoms are conserved from the initial conditions to the equilibrium distribution of species. These constraints are in the form of A_eq x = beq, where x is the vector of mole numbers. CO H2O CO2 H2

Aeq = [ 1    0    1    0   % C balance
        1    1    2    0   % O balance
        0    2    0    2]; % H balance

% equimolar feed of 1 mol H2O and 1 mol CO
beq = [1   % mol C fed
       2   % mol O fed
       2]; % mol H fed

Constraints on bounds

This simply ensures there are no negative mole numbers.

LB = [0 0 0];

Solve the minimization problem

x0 = [0.5 0.5 0.5 0.5]; % initial guesses

options = optimset('Algorithm','sqp');
x = fmincon(@func,x0,[],[],Aeq,beq,LB,[],[],options)

Local minimum found that satisfies the constraints.

Optimization completed because the objective function is non-decreasing in 
feasible directions, to within the default value of the function tolerance,
and constraints are satisfied to within the default value of the constraint tolerance.




x =

    0.4549    0.4549    0.5451    0.5451

Compute mole fractions and partial pressures

The pressures here are in good agreement with the pressures found in Post 1536 . The minor disagreement (in the third or fourth decimal place) is likely due to convergence tolerances in the different algorithms used.

yj = x/sum(x);
Pj = yj*P;

for i=1:numel(x)
    fprintf('%5d%10s%15.4f atm\n',i, species{i}, Pj(i))
end

    1        CO         2.2745 atm
    2       H2O         2.2745 atm
    3       CO2         2.7255 atm
    4        H2         2.7255 atm

Computing equilibrium constants

We can compute the equilibrium constant for the reaction $CO + H_2O \rightleftharpoons CO_2 + H_2$. Compared to the value of K = 1.44 we found at the end of Post 1536 , the agreement is excellent. Note, that to define an equilibrium constant it is necessary to specify a reaction, even though it is not necessary to even consider a reaction to obtain the equilibrium distribution of species!

nuj = [-1 -1 1 1]; % stoichiometric coefficients of the reaction
K = prod(yj.^nuj)

K =

    1.4359


end

% categories: optimization
% tags: Thermodynamics, reaction engineering


Published with MATLAB® 7.12

Read and Post Comments

Finding equilibrium composition by direct minimization of Gibbs free energy on mole numbers

December 25, 2011 at 07:59 PM | categories: optimization | View Comments

Finding equilibrium composition by direct minimization of Gibbs free energy on mole numbers

Finding equilibrium composition by direct minimization of Gibbs free energy on mole numbers

John Kitchin

Adapted from problem 4.5 in Cutlip and Shacham

Contents

Problem setup

Ethane and steam are fed to a steam cracker at a total pressure of 1 atm and at 1000K at a ratio of 4 mol H2O to 1 mol ethane. Estimate the equilibrium distribution of products (CH4, C2H4, C2H2, CO2, CO, O2, H2, H2O, and C2H6).

Solution method: We will construct a Gibbs energy function for the mixture, and obtain the equilibrium composition by minimization of the function subject to elemental mass balance constraints.

function main

R = 0.00198588; % kcal/mol/K
T = 1000; % K

% we store the species names in a cell array since they have different
% lengths (i.e. different number of characters in names).
species = {'CH4' 'C2H4' 'C2H2' 'CO2' 'CO' 'O2' 'H2' 'H2O' 'C2H6'};

%
% $G_^\circ for each species. These are the heats of formation for each
% species.
Gjo = [4.61 28.249 40.604 -94.61 -47.942 0 0 -46.03 26.13]; % kcal/mol

The Gibbs energy of a mixture

We start with $G=\sum\limits_j n_j \mu_j$. Recalling that we define $\mu_j = G_j^\circ + RT \ln a_j$, and in the ideal gas limit, $a_j = y_j P/P^\circ$, and that $y_j = \frac{n_j}{\sum n_j}$. Since in this problem, P = 1 atm, this leads to the function $\frac{G}{RT} = \sum\limits_{j=1}^n n_j\left(\frac{G_j^\circ}{RT} + \ln \frac{n_j}{\sum n_j}\right)$.

    function G = func(nj)
        Enj = sum(nj);
        G = sum(nj.*(Gjo/R/T + log(nj/Enj)));
    end

Linear equality constraints for atomic mass conservation

The total number of each type of atom must be the same as what entered the reactor. These form equality constraints on the equilibrium composition. We express these constraints as: $A_{eq} x = b$ where $x$ is a vector of the moles of each species present in the mixture. CH4 C2H4 C2H2 CO2 CO O2 H2 H2O C2H6

Aeq = [0   0    0    2   1  2  0  1   0       % oxygen balance
       4   4    2    0   0  0  2  2   6       % hydrogen balance
       1   2    2    1   1  0  0  0   2];     % carbon balance

% the incoming feed was 4 mol H2O and 1 mol ethane
beq = [4    % moles of oxygen atoms coming in
       14   % moles of hydrogen atoms coming in
       2];  % moles of carbon atoms coming in

Bounds on mole numbers

No mole number can be negative, so we define a lower bound of zero for each mole number. In principle there are upper bounds too, but it is unecessary to define them here.

LB = [0 0 0 0 0 0 0 0 0]; % no mole numbers less than zero

initial guess for the solver

these are the guesses suggested in the book.

x0 = [1e-3 1e-3 1e-3 0.993 1 1e-4 5.992 1 1e-3]; % initial guess

Setup minimization

fmincon does not solve this problem with the default solver, and warns you to pick another one, either 'interior-point' or 'sqp'. Both work for this problem.

options = optimset('Algorithm','sqp');
[x fval] = fmincon(@func,x0,[],[],Aeq,beq,LB,[],[],options);

for i=1:numel(x)
    fprintf('%5d%10s%10.3g\n',i,species{i},x(i))
end

Local minimum possible. Constraints satisfied.

fmincon stopped because the size of the current step is less than
the default value of the step size tolerance and constraints are 
satisfied to within the default value of the constraint tolerance.



    1       CH4    0.0664
    2      C2H4 8.88e-008
    3      C2H2 2.76e-021
    4       CO2     0.545
    5        CO      1.39
    6        O2 4.56e-021
    7        H2      5.35
    8       H2O      1.52
    9      C2H6  1.6e-007

Verify the equality constraints were met

you can see here the constraints were met.

Aeq*x'
beq

ans =

     4
    14
     2


beq =

     4
    14
     2

Results

Interestingly there is a distribution of products! That is interesting because only steam and ethane enter the reactor, but a small fraction of methane is formed! The main product is hydrogen. The stoichiometry of steam reforming is ideally $C_2H_6 + 4H_2O \rightarrow 2CO_2 + 7 H2$. Even though nearly all the ethane is consumed, we don't get the full yield of hydrogen. It appears that another equilibrium, one between CO, CO2, H2O and H2, may be limiting that, since the rest of the hydrogen is largely in the water. It is also of great importance that we have not said anything about reactions, i.e. how these products were formed. This is in stark contrast to Post 1589 where stoichiometric coefficients were necessary to perform the minimization.

Water gas shift equilibrium

The water gas shift reaction is: $CO + H_2O \rightleftharpoons CO_2 + H_2$. We can compute the Gibbs free energy of the reaction from the heats of formation of each species. Assuming these are the formation energies at 1000K, this is the reaction free energy at 1000K.

G_wgs = Gjo(4) + Gjo(7) - Gjo(5) - Gjo(8)

G_wgs =

   -0.6380

WGS equilibrium constant at 1000K

exp(-G_wgs/R/T)

ans =

    1.3789

Equilibrium constant based on mole numbers

One normally uses activities to define the equilibrium constant. Since there are the same number of moles on each side of the reaction all factors that convert mole numbers to activity, concentration or pressure cancel, so we simply consider the ratio of mole numbers here.

x(4)*x(7)/(x(5)*x(8))

ans =

    1.3789

Clearly, there is an equilibrium between these species that prevents the complete reaction of steam reforming.

Summary

This is an appealing way to minimize the Gibbs energy of a mixture. No assumptions about reactions are necessary, and the constraints are easy to identify. The Gibbs energy function is especially easy to code.

end

% categories: optimization
% tags: Thermodynamics, reaction engineering
% post_id = 1630; %delete this line to force new post;
% permaLink = http://matlab.cheme.cmu.edu/2011/12/25/finding-equilibrium-composition-by-direct-minimization-of-gibbs-free-energy-on-mole-numbers/;


Published with MATLAB® 7.13

Read and Post Comments

The Gibbs free energy of a reacting mixture and the equilibrium composition

December 20, 2011 at 05:25 PM | categories: optimization | View Comments

The Gibbs free energy of a reacting mixture and the equilibrium composition

The Gibbs free energy of a reacting mixture and the equilibrium composition

John Kitchin

Adapted from Chapter 3 in Rawlings and Ekerdt.

Contents

Introduction

In this post we derive the equations needed to find the equilibrium composition of a reacting mixture. We use the method of direct minimization of the Gibbs free energy of the reacting mixture.

The Gibbs free energy of a mixture is defined as $G = \sum\limits_j \mu_j n_j$ where $\mu_j$ is the chemical potential of species $j$, and it is temperature and pressure dependent, and $n_j$ is the number of moles of species $j$.

We define the chemical potential as $\mu_j = G_j^\circ + RT\ln a_j$, where $G_j^\circ$ is the Gibbs energy in a standard state, and $a_j$ is the activity of species $j$ if the pressure and temperature are not at standard state conditions.

If a reaction is occurring, then the number of moles of each species are related to each other through the reaction extent $\epsilon$ and stoichiometric coefficients: $n_j = n_{j0} + \nu_j \epsilon$. Note that the reaction extent has units of moles.

Combining these three equations and expanding the terms leads to:

$$G = \sum\limits_j n_{j0}G_j^\circ +\sum\limits_j \nu_j G_j^\circ \epsilon +RT\sum\limits_j(n_{j0} + \nu_j\epsilon)\ln a_j $$

The first term is simply the initial Gibbs free energy that is present before any reaction begins, and it is a constant. It is difficult to evaluate, so we will move it to the left side of the equation in the next step, because it does not matter what its value is since it is a constant. The second term is related to the Gibbs free energy of reaction: $\Delta_rG = \sum\limits_j \nu_j G_j^\circ$. With these observations we rewrite the equation as:

$$G - \sum\limits_j n_{j0}G_j^\circ = \Delta_rG \epsilon +RT\sum\limits_j(n_{j0} + \nu_j\epsilon)\ln a_j $$

Now, we have an equation that allows us to compute the change in Gibbs free energy as a function of the reaction extent, initial number of moles of each species, and the activities of each species. This difference in Gibbs free energy has no natural scale, and depends on the size of the system, i.e. on $n_{j0}$. It is desirable to avoid this, so we now rescale the equation by the total initial moles present, $n_{T0}$ and define a new variable $\epsilon' = \epsilon/n_{T0}$, which is dimensionless. This leads to:

$$ \frac{G - \sum\limits_j n_{j0}G_j^\circ}{n_{T0}} = \Delta_rG \epsilon' + RT \sum\limits_j(y_{j0} + \nu_j\epsilon')\ln a_j $$

where $y_{j0}$ is the initial mole fraction of species $j$ present. The mole fractions are intensive properties that do not depend on the system size. Finally, we need to address $a_j$. For an ideal gas, we know that $A_j = \frac{y_j P}{P^\circ}$, where the numerator is the partial pressure of species $j$ computed from the mole fraction of species $j$ times the total pressure. To get the mole fraction we note:

$$y_j = \frac{n_j}{n_T} = \frac{n_{j0} + \nu_j \epsilon}{n_{T0} + \epsilon \sum\limits_j \nu_j} = \frac{y_{j0} + \nu_j \epsilon'}{1 + \epsilon'\sum\limits_j \nu_j} $$

This finally leads us to an equation that we can evaluate as a function of reaction extent:

$$ \frac{G - \sum\limits_j n_{j0}G_j^\circ}{n_{T0}} = \widetilde{\widetilde{G}} = \Delta_rG \epsilon' + RT\sum\limits_j(y_{j0} + \nu_j\epsilon') \ln\left(\frac{y_{j0}+\nu_j\epsilon'}{1+\epsilon'\sum\limits_j\nu_j} \frac{P}{P^\circ}\right) $$

we use a double tilde notation to distinguish this quantity from the quantity derived by Rawlings and Ekerdt which is further normalized by a factor of $RT$. This additional scaling makes the quantities dimensionless, and makes the quantity have a magnitude of order unity, but otherwise has no effect on the shape of the graph.

Finally, if we know the initial mole fractions, the initial total pressure, the Gibbs energy of reaction, and the stoichiometric coefficients, we can plot the scaled reacting mixture energy as a function of reaction extent. At equilibrium, this energy will be a minimum. We consider the example in Rawlings and Ekerdt where isobutane (I) reacts with 1-butene (B) to form 2,2,3-trimethylpentane (P). The reaction occurs at a total pressure of 2.5 atm at 400K, with equal molar amounts of I and B. The standard Gibbs free energy of reaction at 400K is -3.72 kcal/mol. Compute the equilibrium composition.

function main

clear all; close all

u = cmu.units;
R = 8.314*u.J/u.mol/u.K;
P = 2.5*u.atm;
Po = 1*u.atm; % standard state pressure
T = 400*u.K;

Grxn = -3.72*u.kcal/u.mol; %reaction free energy at 400K
yi0 = 0.5; yb0 = 0.5; yp0 = 0.0; % initial mole fractions

yj0 = [yi0 yb0 yp0];
nu_j = [-1 -1 1]; % stoichiometric coefficients

Define a function for $\frac{G - \sum\limits_j n_{j0}G_j^\circ}{n_{T0}}$

We will call this difference Gwigglewiggle.

    function diffG = Gfunc(extentp)
        diffG = Grxn*extentp;

        sum_nu_j = sum(nu_j);
        for i = 1:length(yj0)
            x1 = yj0(i) + nu_j(i)*extentp;
            x2 = x1./(1+extentp*sum_nu_j);
            diffG = diffG + R*T*x1.*log(x2*P/Po);
        end
    end

Plot function

There are bounds on how large $\epsilon'$ can be. Recall that $n_j = n_{j0} + \nu_j \epsilon$, and that $n_j \ge 0$. Thus, $\epsilon_{max} = -n_{j0}/\nu_j$, and the maximum value that $\epsilon'$ can have is therefore $-y_{j0}/\nu_j$ where $y_{j0}>0$. When there are multiple species, you need the smallest $epsilon'_{max}$ to avoid getting negative mole numbers.

epsilonp_max = min(-yj0(yj0>0)./nu_j(yj0>0))
epsilonp = linspace(0,epsilonp_max,1000);
plot(epsilonp,Gfunc(epsilonp))
xlabel('\epsilon''')
ylabel('Gwigglewiggle')

epsilonp_max =

    0.5000

Find extent that minimizes Gwigglewiggle

fminbnd is not defined for the units class, so we wrap the function in a function handle that makes it output a double. Alternatively, we could make the function fully dimensionless. We save that for a later post.

f = @(x) double(Gfunc(x))
epsilonp_eq = fminbnd(f,0,epsilonp_max)

hold on
plot(epsilonp_eq,Gfunc(epsilonp_eq),'ro')

f = 

    @(x)double(Gfunc(x))


epsilonp_eq =

    0.4696

Compute equilibrium mole fractions

yi = (yi0 + nu_j(1)*epsilonp_eq)/(1 + epsilonp_eq*sum(nu_j))
yb = (yb0 + nu_j(2)*epsilonp_eq)/(1 + epsilonp_eq*sum(nu_j))
yp = (yp0 + nu_j(3)*epsilonp_eq)/(1 + epsilonp_eq*sum(nu_j))

yi =

    0.0572


yb =

    0.0572


yp =

    0.8855

alternative way to compute mole fractions

y_j = (yj0 + nu_j*epsilonp_eq)/(1 + epsilonp_eq*sum(nu_j))

y_j =

    0.0572    0.0572    0.8855

Check for consistency with the equilibrium constant

K = exp(-Grxn/R/T)

K =

  108.1300

$K = \frac{a_P}{a_I a_B} = \frac{y_p P/P^\circ}{y_i P/P^\circ y_b P/P^\circ} = \frac{y_P}{y_i y_b}\frac{P^\circ}{P}$.

yp/(yi*yb)*Po/P

ans =

  108.0702

an alternative form of the equilibrium constant

We can express the equilibrium constant like this :$K = \prod\limits_j a_j^{\nu_j}$, and compute it with a single command.

prod((y_j*P/Po).^nu_j)

ans =

  108.0702

These results are very close, and only disagree because of the default tolerance used in identifying the minimum of our function. you could tighten the tolerances by setting options to the fminbnd function.

Summary

In this post we derived an equation for the Gibbs free energy of a reacting mixture and used it to find the equilibrium composition. In future posts we will examine some alternate forms of the equations that may be more useful in some circumstances.

end

% categories: optimization
% tags: thermodynamics, reaction engineering


Published with MATLAB® 7.13

Read and Post Comments

Conservation of mass in chemical reactions

December 19, 2011 at 02:03 AM | categories: uncategorized | View Comments

Conservation of mass in chemical reactions

Conservation of mass in chemical reactions

John Kitchin

Atoms cannot be destroyed in non-nuclear chemical reactions, hence it follows that the same number of atoms entering a reactor must also leave the reactor. The atoms may leave the reactor in a different molecular configuration due to the reaction, but the total mass leaving the reactor must be the same. Here we look at a few ways to show this.

Contents

Water gas shift reaction

We consider the water gas shift reaction : % $CO + H_2O \rightleftharpoons H_2 + CO_2$. We can illustrate the conservation of mass with the following equation: $\bf{\nu}\bf{M}=\bf{0}$. Where $\bf{\nu}$ is the stoichiometric coefficient vector and $\bf{M}$ is a column vector of molecular weights. For simplicity, we use pure isotope molecular weights, and not the isotope-weighted molecular weights.

nu = [-1 -1 1 1];
M = [28; 18; 2; 44];
nu*M

ans =

     0

Atomic mass balances

For any balanced chemical equation, there are the same number of each kind of atom on each side of the equation. Since the mass of each atom is unchanged with reaction, that means the mass of all the species that are reactants must equal the mass of all the species that are products! Here we look at the n C O H

reactants = [-1 -2 -2]
products  = [ 1  2  2]
M = [12.011; 15.999; 1.0079]

reactants =

    -1    -2    -2


products =

     1     2     2


M =

   12.0110
   15.9990
    1.0079

Now if we add the mass of reactants and products, it should sum to zero (since we used the negative sign for stoichiometric coefficients of reactants).

products*M + reactants*M

ans =

     0

Summary

That's all there is to mass conservation with reactions. Nothing changes if there are lots of reactions, as long as each reaction is properly balanced, and none of them are nuclear reactions!

tags: reaction engineering


Published with MATLAB® 7.12

Read and Post Comments

Water gas shift equilibria via the NIST Webbook

December 13, 2011 at 01:49 AM | categories: miscellaneous, nonlinear algebra | View Comments

Water gas shift equilibria via the NIST Webbook

Water gas shift equilibria via the NIST Webbook

John Kitchin

The NIST webbook provides parameterized models of the enthalpy, entropy and heat capacity of many molecules. In this example, we will examine how to use these to compute the equilibrium constant for the water gas shift reaction $CO + H_2O \rightleftharpoons CO_2 + H_2$ in the temperature range of 500K to 1000K.

Parameters are provided for:

Contents

   Cp = heat capacity (J/mol*K)
   H° = standard enthalpy (kJ/mol)
   S° = standard entropy (J/mol*K)

with models in the form: $Cp^\circ = A + B*t + C*t^2 + D*t^3 + E/t^2$

$H^\circ - H^\circ_{298.15}= A*t + B*t^2/2 + C*t^3/3 + D*t^4/4 - E/t + F - H$

$S^\circ = A*ln(t) + B*t + C*t^2/2 + D*t^3/3 - E/(2*t^2) + G$

where $t=T/1000$, and $T$ is the temperature in Kelvin. We can use this data to calculate equilibrium constants in the following manner. First, we have heats of formation at standard state for each compound; for elements, these are zero by definition, and for non-elements, they have values available from the NIST webbook. There are also values for the absolute entropy at standard state. Then, we have an expression for the change in enthalpy from standard state as defined above, as well as the absolute entropy. From these we can derive the reaction enthalpy, free energy and entropy at standard state, as well as at other temperatures.

We will examine the water gas shift enthalpy, free energy and equilibrium constant from 500K to 1000K, and finally compute the equilibrium composition of a gas feed containing 5 atm of CO and H2 at 1000K.

close all; clear all; clc;

T = linspace(500,1000); % degrees K
t = T/1000;

Setup equations for each species

First we enter in the parameters and compute the enthalpy and entropy for each species.

Hydrogen

link

% T = 298-1000K valid temperature range
A =  33.066178;
B = -11.363417;
C =  11.432816;
D = -2.772874;
E = -0.158558;
F = -9.980797;
G =  172.707974;
H =  0.0;

Hf_29815_H2 = 0.0; % kJ/mol
S_29815_H2 = 130.68; % J/mol/K

dH_H2 = A*t + B*t.^2/2 + C*t.^3/3 + D*t.^4/4 - E./t + F - H;
S_H2 = (A*log(t) + B*t + C*t.^2/2 + D*t.^3/3 - E./(2*t.^2) + G);

H2O

link Note these parameters limit the temperature range we can examine, as these parameters are not valid below 500K. There is another set of parameters for lower temperatures, but we do not consider them here.

% 500-1700 K valid temperature range
A =   30.09200;
B =   6.832514;
C =   6.793435;
D =  -2.534480;
E =   0.082139;
F =  -250.8810;
G =   223.3967;
H =  -241.8264;

Hf_29815_H2O = -241.83; %this is Hf.
S_29815_H2O = 188.84;

dH_H2O = A*t + B*t.^2/2 + C*t.^3/3 + D*t.^4/4 - E./t + F - H;
S_H2O = (A*log(t) + B*t + C*t.^2/2 + D*t.^3/3 - E./(2*t.^2) + G);

CO

link

% 298. - 1300K valid temperature range
A =   25.56759;
B =   6.096130;
C =   4.054656;
D =  -2.671301;
E =   0.131021;
F =  -118.0089;
G =   227.3665;
H = -110.5271;

Hf_29815_CO = -110.53; %this is Hf kJ/mol.
S_29815_CO = 197.66;

dH_CO = A*t + B*t.^2/2 + C*t.^3/3 + D*t.^4/4 - E./t + F - H;
S_CO = (A*log(t) + B*t + C*t.^2/2 + D*t.^3/3 - E./(2*t.^2) + G);

CO2

link

% 298. - 1200.K valid temperature range
A =   24.99735;
B =   55.18696;
C =  -33.69137;
D =   7.948387;
E =  -0.136638;
F =  -403.6075;
G =   228.2431;
H =  -393.5224;

Hf_29815_CO2 = -393.51; %this is Hf.
S_29815_CO2 = 213.79;

dH_CO2 = A*t + B*t.^2/2 + C*t.^3/3 + D*t.^4/4 - E./t + F - H;
S_CO2 = (A*log(t) + B*t + C*t.^2/2 + D*t.^3/3 - E./(2*t.^2) + G);

Standard state heat of reaction

We compute the enthalpy and free energy of reaction at 298.15 K for the following reaction $CO + H2O \rightleftharpoons H2 + CO2$.

Hrxn_29815 = Hf_29815_CO2 + Hf_29815_H2 - Hf_29815_CO - Hf_29815_H2O;
Srxn_29815 = S_29815_CO2 + S_29815_H2 - S_29815_CO - S_29815_H2O;
Grxn_29815 = Hrxn_29815 - 298.15*(Srxn_29815)/1000;

sprintf('deltaH = %1.2f',Hrxn_29815)
sprintf('deltaG = %1.2f',Grxn_29815)

ans =

deltaH = -41.15


ans =

deltaG = -28.62

Correcting $\Delta H$ and $\Delta G$

we have to correct for temperature change away from standard state. We only correct the enthalpy for this temperature change. The correction looks like this:

$$ \Delta H_{rxn}(T) = \Delta H_{rxn}(T_{ref}) + \sum_i \nu_i (H_i(T)-H_i(T_{ref}))$$

Where $\nu_i$ are the stoichiometric coefficients of each species, with appropriate sign for reactants and products, and $(H_i(T)-H_i(T_{ref})$ is precisely what is calculated for each species with the equations

Hrxn = Hrxn_29815 + dH_CO2 + dH_H2 - dH_CO - dH_H2O;

The entropy is on an absolute scale, so we directly calculate entropy at each temperature. Recall that H is in kJ/mol and S is in J/mol/K, so we divide S by 1000 to make the units match.

Grxn = Hrxn - T.*(S_CO2 + S_H2 - S_CO - S_H2O)/1000;

Plot how the $\Delta G$ varies with temperature

over this temperature range the reaction is exothermic, although near 1000K it is just barely exothermic. At higher temperatures we expect the reaction to become endothermic.

figure; hold all
plot(T,Grxn)
plot(T,Hrxn)
xlabel('Temperature (K)')
ylabel('(kJ/mol)')
legend('\Delta G_{rxn}', '\Delta H_{rxn}', 'location','best')

Equilibrium constant calculation

Note the equilibrium constant starts out high, i.e. strongly favoring the formation of products, but drops very quicky with increasing temperature.

R = 8.314e-3; %kJ/mol/K
K = exp(-Grxn/R./T);

figure
plot(T,K)
xlim([500 1000])
xlabel('Temperature (K)')
ylabel('Equilibrium constant')

Equilibrium yield of WGS

Now let's suppose we have a reactor with a feed of H2O and CO at 10atm at 1000K. What is the equilibrium yield of H2? Let $\epsilon$ be the extent of reaction, so that $F_i = F_{i,0} + \nu_i \epsilon$. For reactants, $\nu_i$ is negative, and for products, $\nu_i$ is positive. We have to solve for the extent of reaction that satisfies the equilibrium condition.

A = CO
B = H2O
C = H2
D = CO2
Pa0 = 5; Pb0 = 5; Pc0 = 0; Pd0 = 0;  % pressure in atm
R = 0.082;
Temperature = 1000;

% we can estimate the equilibrium like this. We could also calculate it
% using the equations above, but we would have to evaluate each term. Above
% we simple computed a vector of enthalpies, entropies, etc...
K_Temperature = interp1(T,K,Temperature);

% If we let X be fractional conversion then we have $C_A = C_{A0}(1-X)$,
% $C_B = C_{B0}-C_{A0}X$, $C_C = C_{C0}+C_{A0}X$, and $C_D =
% C_{D0}+C_{A0}X$. We also have $K(T) = (C_C C_D)/(C_A C_B)$, which finally
% reduces to $0 = K(T) - Xeq^2/(1-Xeq)^2$ under these conditions.

f = @(X) K_Temperature - X^2/(1-X)^2;

Xeq = fzero(f,[1e-3 0.999]);
sprintf('The equilibrium conversion for these feed conditions is: %1.2f',Xeq)

ans =

The equilibrium conversion for these feed conditions is: 0.55

Compute gas phase pressures of each species

Since there is no change in moles for this reaction, we can directly calculation the pressures from the equilibrium conversion and the initial pressure of gases. you can see there is a slightly higher pressure of H2 and CO2 than the reactants, consistent with the equilibrium constant of about 1.44 at 1000K. At a lower temperature there would be a much higher yield of the products. For example, at 550K the equilibrium constant is about 58, and the pressure of H2 is 4.4 atm due to a much higher equilibrium conversion of 0.88.

P_CO = Pa0*(1-Xeq)
P_H2O = Pa0*(1-Xeq)
P_H2 = Pa0*Xeq
P_CO2 = Pa0*Xeq

P_CO =

    2.2748


P_H2O =

    2.2748


P_H2 =

    2.7252


P_CO2 =

    2.7252

Compare the equilibrium constants

We can compare the equilibrium constant from the Gibbs free energy and the one from the ratio of pressures. They should be the same!

K_Temperature
(P_CO2*P_H2)/(P_CO*P_H2O)

K_Temperature =

    1.4352


ans =

    1.4352

Summary

The NIST Webbook provides a plethora of data for computing thermodynamic properties. It is a little tedious to enter it all into Matlab, and a little tricky to use the data to estimate temperature dependent reaction energies. A limitation of the Webbook is that it does not tell you have the thermodynamic properties change with pressure. Luckily, those changes tend to be small.

% categories: Miscellaneous, nonlinear algebra
% tags: thermodynamics, reaction engineering


Published with MATLAB® 7.12

Read and Post Comments
Next Page »