The Kitchin Research Group

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

yj = N / np.sum(N)
Pj = yj * P

for s, y, p in zip(species, yj, Pj):
    print '{0:10s}: {1:1.2f} {2:1.2f}'.format(s, y, p)

>>> >>> ... ... CO        : 0.23 2.28
H2O       : 0.23 2.28
CO2       : 0.27 2.72
H2        : 0.27 2.72

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 = np.array([-1, -1, 1, 1])  # stoichiometric coefficients of the reaction
K = np.prod(yj**nuj)
print K

>>> 1.43446295961

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)\).

import numpy as np

def func(nj):
    nj = np.array(nj)
    Enj = np.sum(nj);
    G = np.sum(nj * (Gjo / R / T + np.log(nj / Enj)))
    return G

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} n = b\) where \(n\) is a vector of the moles of each species present in the mixture. CH4 C2H4 C2H2 CO2 CO O2 H2 H2O C2H6

Aeq = np.array([[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 = np.array([4,  # moles of oxygen atoms coming in
                14, # moles of hydrogen atoms coming in
                2]) # moles of carbon atoms coming in

def ec1(n):
    'equality constraint'
    return np.dot(Aeq, n) - beq

def ic1(n):
    '''inequality constraint
       all n>=0
    '''   
    return n

Now we solve the problem.

# initial guess suggested in the example
n0 = [1e-3, 1e-3, 1e-3, 0.993, 1.0, 1e-4, 5.992, 1.0, 1e-3] 

n0 = [0.066, 8.7e-08, 2.1e-14, 0.545, 1.39, 5.7e-14, 5.346, 1.521, 1.58e-7]

from scipy.optimize import fmin_slsqp

X = fmin_slsqp(func, n0, f_eqcons=ec1,f_ieqcons=ic1, iter=300, acc=1e-12)

for s,x in zip(species, X):
    print '{0:10s} {1:1.4g}'.format(s, x)

# check that constraints were met
print np.dot(Aeq, X) - beq
print np.all( np.abs( np.dot(Aeq, X) - beq) < 1e-12)

>>> >>> >>> >>> >>> >>> Optimization terminated successfully.    (Exit mode 0)
            Current function value: -104.403951524
            Iterations: 16
            Function evaluations: 193
            Gradient evaluations: 15
>>> ... ... CH4        0.06644
C2H4       9.48e-08
C2H2       1.487e-13
CO2        0.545
CO         1.389
O2         3.096e-13
H2         5.346
H2O        1.521
C2H6       1.581e-07
... [  0.00000000e+00   0.00000000e+00   4.44089210e-16]
True

I found it necessary to tighten the accuracy parameter to get pretty good matches to the solutions found in Matlab. It was also necessary to increase the number of iterations. Even still, not all of the numbers match well, especially the very small numbers. You can, however, see that the constraints were satisfied pretty well.

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 do not 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.

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[3] + Gjo[6] - Gjo[4] - Gjo[7]
print G_wgs

K = np.exp(-G_wgs / (R*T))
print K

-0.638
>>> >>> 1.37887528109

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.

print (X[3] * X[6]) / (X[4] * X[7])

1.37887525547

This is very close to the equilibrium constant computed above.

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

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.

We have saturated liquid here, and we get the thermodynamic properties for the given temperature.

from iapws import IAPWS97

T1 = 100 + 273.15 #in K

sat_liquid1  = IAPWS97(T=T1, x=0) # x is the steam quality. 0 = liquid

P1 = sat_liquid1.P
s1 = sat_liquid1.s
h1 = sat_liquid1.h
v1 = sat_liquid1.v

The final pressure is given, and we need to compute the new temperatures, and enthalpy.

P2 = 3.0 # MPa
s2 = s1 # this is what isentropic means

sat_liquid2 = IAPWS97(P=P2, s=s1)
T2, = sat_liquid2.T
h2 = sat_liquid2.h

# work done to compress liquid. This is an approximation, since the
# volume does change a little with pressure, but the overall work here
# is pretty small so we neglect the volume change.
WdotP = v1*(P2 - P1);
print
print('The compressor work is: {0:1.4f} kJ/kg'.format(WdotP))

>>> >>> >>> >>> >>> >>> ... ... ... >>>
The compressor work is: 0.0030 kJ/kg

The compression work is almost negligible. This number is 1000 times smaller than we computed with Xsteam. I wonder what the units of v1 actually are.

T3 = 600 + 273.15 # K
P3 = P2 # definition of isobaric
steam = IAPWS97(P=P3, T=T3)

h3 = steam.h
s3 = steam.s

Qb, = h3 - h2 # heat required to make the steam

print
print 'The boiler heat duty is: {0:1.2f} kJ/kg'.format(Qb)

>>> >>> >>> >>> >>> >>> >>> >>>
The boiler heat duty is: 3260.69 kJ/kg

steam =  IAPWS97(P=P1, s=s3)
T4, = steam.T
h4 = steam.h
s4 = s3 # isentropic
Qc, = h4 - h1 # work required to cool from T4 to T1
print 
print 'The condenser heat duty is {0:1.2f} kJ/kg'.format(Qc)

>>> >>> >>> >>>
The condenser heat duty is 2317.00 kJ/kg

WdotTurbine, = h4 - h3 # work extracted from the expansion
print('The turbine work is: {0:1.2f} kJ/kg'.format(WdotTurbine))

The turbine work is: -946.71 kJ/kg

This is a ratio of the work put in to make the steam, and the net work obtained from the turbine. The answer here agrees with the efficiency calculated in Sandler on page 135.

eta = -(WdotTurbine - WdotP) / Qb
print('The overall efficiency is {0:1.2%}.'.format(eta))

The overall efficiency is 29.03%.

The IAPWS module makes it pretty easy to generate figures of the steam tables. Here we generate an entropy-Temperature graph. We do this to illustrate the path of the Rankine cycle. We need to compute the values of steam entropy for a range of pressures and temperatures.

import numpy as np
import matplotlib.pyplot as plt

plt.figure()
plt.clf()
T = np.linspace(300, 372+273, 200) # range of temperatures
for P in [0.1, 1, 2, 5, 10, 20]: #MPa
    steam = [IAPWS97(T=t, P=P) for t in T]
    S = [s.s for s in steam]
    plt.plot(S, T, 'k-')

# saturated vapor and liquid entropy lines
svap = [s.s for s in [IAPWS97(T=t, x=1) for t in T]]
sliq = [s.s for s in [IAPWS97(T=t, x=0) for t in T]]

plt.plot(svap, T, 'r-')
plt.plot(sliq, T, 'b-')

plt.xlabel('Entropy (kJ/(kg K)')
plt.ylabel('Temperature (K)')
plt.savefig('images/iawps-steam.png')

>>> >>> <matplotlib.figure.Figure object at 0x000000000638BC18>
>>> >>> ... ... ... ... [<matplotlib.lines.Line2D object at 0x0000000007F9C208>]
[<matplotlib.lines.Line2D object at 0x0000000007F9C400>]
[<matplotlib.lines.Line2D object at 0x0000000007F9C8D0>]
[<matplotlib.lines.Line2D object at 0x0000000007F9CD30>]
[<matplotlib.lines.Line2D object at 0x0000000007F9E1D0>]
[<matplotlib.lines.Line2D object at 0x0000000007F9E630>]
... >>> >>> >>> [<matplotlib.lines.Line2D object at 0x0000000001FDCEB8>]
[<matplotlib.lines.Line2D object at 0x0000000007F9EA90>]
>>> <matplotlib.text.Text object at 0x0000000007F7BE48>
<matplotlib.text.Text object at 0x0000000007F855F8>

We can plot our Rankine cycle path like this. We compute the entropies along the non-isentropic paths.

T23 = np.linspace(T2, T3)
S23 = [s.s for s in [IAPWS97(P=P2, T=t) for t in T23]]

T41 = np.linspace(T4, T1 - 0.01) # subtract a tiny bit to make sure we get a liquid
S41 = [s.s for s in [IAPWS97(P=P1, T=t) for t in T41]]

And then we plot the paths.

plt.plot([s1, s2], [T1, T2], 'r-', lw=4) # Path 1 to 2
plt.plot(S23, T23, 'b-', lw=4) # path from 2 to 3 is isobaric
plt.plot([s3, s4], [T3, T4], 'g-', lw=4) # path from 3 to 4 is isentropic
plt.plot(S41, T41, 'k-', lw=4) # and from 4 to 1 is isobaric
plt.savefig('images/iawps-steam-2.png')
plt.savefig('images/iawps-steam-2.svg')

[<matplotlib.lines.Line2D object at 0x0000000008350908>]
[<matplotlib.lines.Line2D object at 0x00000000083358D0>]
[<matplotlib.lines.Line2D object at 0x000000000835BEB8>]
[<matplotlib.lines.Line2D object at 0x0000000008357160>]

This was an interesting exercise. On one hand, the tedium of interpolating the steam tables is gone. On the other hand, you still have to know exactly what to ask for to get an answer that is correct. The iapws interface is a little clunky, and takes some getting used to. It does not seem as robust as the Xsteam module I used in Matlab.