J.C. Kantor ([email protected])

In [2]:

```
#Initializations
from IPython.core.display import HTML
HTML(open("../styles/custom.css", "r").read())
```

Out[2]:

The steam reforming of methane

$$ CH_4(g) + H_2O(g) \longrightarrow CO(g) + 3\,H_2(g) $$is an endothermic reaction for the production of hydrogen. The reaction is normally conducted at relatively high temperatures and pressures over a heterogeneous catalyst. Assuming the reaction is carried out with 20% excess water at 20 atm, solve for the equilibrium extent of reaction as a function of temperature. Then plot methane predicted conversion as a function of operating temperature.

For an ideal mixture of gases at equilibrium,

$$K_{rxn}(T) = \prod_{c=1}^C \left(y_cP\right)^{\nu_c}$$or, after taking logarithms,

$$\ln K_{rxn}(T) = \sum_{c=1}^C \nu_c \ln\left(y_cP\right)$$where $K_{rxn}(T)$ is the equilibrium constant, $y_c$ is the mole fraction of component $c$, $\nu_c$ is the stoichiometric coefficient of component $c$, $P$ and $T$ are pressure and temperature, respectively.

The strategy for solving this problem is to use the Van't Hoff equation to write a Python function to compute the left-hand side of the equation as a function of temperature, and the component material balances to solve for the quotient on the right-hand side as a function of the extent of reaction. Then, given a temperature and pressure, the equation is solved for the equilibrium extent of reaction using a root-finding algorithm.

In [1]:

```
pylab.rcParams['figure.figsize'] = (10,4)
pylab.rcParams['font.size'] = 12
pylab.rcParams['lines.linewidth'] = 1.6
```

In [2]:

```
# Reaction Conditions
P = 20
Tmin = 800
Tmax = 1200
# Reactor Feed Specifications
nIn = dict()
nIn['CH4'] = 1
nIn['H2O'] = 1.2*nIn['CH4']
nIn['CO'] = 0
nIn['H2'] = 0
# List of Components
C = nIn.keys()
```

We begin by gathering the molar enthalpy and Gibb's free energy of the species participating in the reaction at standard conditions.

Species | $\Delta\hat{G}_f^{\circ}$ [kJ/gmol] |
$\Delta\hat{H}_f^{\circ}$ [kJ/gmol] |
---|---|---|

Carbon Monoxide (g) | -137.27 | -110.53 |

Hydrogen (g) | 0 | 0 |

Methane (g) | -50.49 | -74.52 |

Water (g) | -228.59 | -241.83 |

In [3]:

```
# Enthalpy and Free Energy of Formation
Gf = dict()
Hf = dict()
Gf['CO'] = -137270
Hf['CO'] = -110530
Gf['H2'] = 0
Hf['H2'] = 0
Gf['CH4'] = -50490
Hf['CH4'] = -74520
Gf['H2O'] = -228590
Hf['H2O'] = -241830
```

The enthalpy and Gibb's free energy per mole extent of reaction are given by

$$\begin{align*} \Delta\hat{H}_{rxn}^{\circ} & = \sum_{c=1}^C \nu_c \Delta\hat{H}_{f,c}^{\circ} \\ \Delta\hat{G}_{rxn}^{\circ} & = \sum_{c=1}^C \nu_c \Delta\hat{G}_{f,c}^{\circ} \end{align*}$$where $\nu_c$ is the stoichiometric coefficient for species $c$.

In [4]:

```
# Reaction Stoichiometry
nu = dict()
nu['CO'] = +1
nu['H2'] = +3
nu['CH4'] = -1
nu['H2O'] = -1
Hr = 0.0
for c in C:
Hr += nu[c]*Hf[c]
print "Hr = {:7.3f} kJ/gmol".format(Hr/1000)
Gr = 0.0
for c in C:
Gr += nu[c]*Gf[c]
print "Gr = {:7.3f} kJ/gmol".format(Gr/1000)
```

The equilibrium constant as a function of temperature can be estimated using the Van't Hoff equation.

$$ \ln K_r(T) = -\frac{1}{R}\left[\frac{\Delta\hat{G}_{rxn}^{\circ} - \Delta\hat{H}_{rxn}^{\circ}}{298} + \frac{\Delta\hat{H}_{rxn}^{\circ}}{T}\right] $$This is implemented in Python as a lambda function of $T$

In [5]:

```
R = 8.314 # J/K/gmol
lnKr = lambda T: -((Gr-Hr)/298 + Hr/T)/R
```

which can be plotted over the temperature range of interest.

In [6]:

```
T = linspace(Tmin,Tmax,100)
plot(T,map(lnKr,T))
xlabel('Temperature [K]')
ylabel('$\ln(K_{r})$')
title('Van\'t Hoff Equation')
grid()
```

Given the problem data, the material balances can be solved to determine the molar flows at the outlet of the reactor as functions of the molar extent of reaction $\dot{\xi}$

$$\begin{align*} \dot{n}_{out,CO} & = \dot{n}_{in,CO} + \nu_{CO}\dot{\xi} \\ \dot{n}_{out,H_2} & = \dot{n}_{in,H_2} + \nu_{H_2}\dot{\xi} \\ \dot{n}_{out,CH_4} & = \dot{n}_{in,CH_4} + \nu_{CH_4}\dot{\xi} \\ \dot{n}_{out,H_2O} & = \dot{n}_{in,H_2O} + \nu_{H_2O}\dot{\xi} \end{align*}$$In [7]:

```
nOut = dict()
for c in C:
nOut[c] = lambda x,c=c: nIn[c] + nu[c]*x
```

In [8]:

```
xMax = Inf
for c in C:
if nu[c] < 0:
xMax = min(xMax,-nIn[c]/nu[c])
print 'Maximum molar extent of reaction = {:.2f}'.format(xMax)
x = linspace(0,xMax)
for c in C:
plot(x,map(nOut[c],x))
legend(C,'best')
xlabel('Molar Extent of Reaction')
ylabel('$n_{Out}$')
title('Molar flows out the Reactor as Function of the Extent of Reaction')
grid()
```

The composition of the reactor outlet gases is given by

$$ y_{n}(\dot{\xi}) = \frac{\dot{n}_{out,n}(\dot{\xi})}{\dot{n}_{Total}(\dot{\xi})} $$where the total molar flow is sum of the component molar flows.

$$\begin{align*} \dot{n}_{Total} & = \sum_{c=1}^C \dot{n}_{out,c} \\ & = \sum_{c=1}^C \dot{n}_{in,c} + \left(\sum_{c=1}^C \nu_c\right)\dot{\xi} \end{align*}$$The total molar flow is implement as a Python function.

In [9]:

```
def nTotal(x):
nTotal = 0
for c in C:
nTotal += nOut[c](x)
return nTotal
```

In [10]:

```
y = dict()
for c in C:
y[c] = lambda x,c=c: nOut[c](x)/nTotal(x)
```

Plotting composition as a function of extent of reaction.

In [11]:

```
x = linspace(0,1)
for c in C:
plot(x,map(y[c],x))
legend(C,'best')
xlabel('Molar Extent of Reaction')
ylabel('$y_c$')
title('Composition as a Function of Molar Extent of Reaction')
grid()
```

For a gas phase reaction involving a mixture of ideal gases, we define a reaction quotient as

$$K_a(\dot{\xi}) = \prod_{c=1}^C \left(y_c(\dot{\xi})P\right)^{\nu_c}$$where $y_c(\dot{\xi})P$ is the partial pressure of component $c$. Taking the logarithm

$$ \ln K_a(\dot{\xi}) = \sum_{c=1}^C \nu_c\ln(y_c(\dot{\xi})P)$$This is implemented as a Python function.

In [12]:

```
def lnKa(x):
lnKa = 0;
for c in C:
lnKa += nu[c]*log(P*y[c](x))
return lnKa
```

In [13]:

```
x = linspace(0.001,1)
plot(x,map(lnKa,x))
xlabel('Molar Extent of Reaction')
ylabel('$\ln(Ka)$')
title('Reaction Activity Quotient');
grid()
```

For a given temperature and pressure, the solution for the equilibrium extent of reaction is the value of $\dot{\xi}$ for which $K_r(T) = K_a(\dot{\xi})$ which can be written in terms of logarithms as

$$\ln K_r(T) = \ln K_a(\dot{\xi})$$Plotting these functions side-by-side shows a simple graphical technique for finding solutions for $\dot{\xi}$.

In [14]:

```
figure(figsize=(12,5))
subplot(1,2,1)
plot(T,map(lnKr,T))
ylim(-5,10)
xlabel('Temperature [K]')
ylabel('$\ln(Kr)$')
title('Equibrium Constant')
grid()
subplot(1,2,2)
plot(x,map(lnKa,x))
ylim(-5,10)
xlabel('Molar Extent of Reaction')
ylabel('$ln(Ka)$')
title('Equilibrium Activity Quotient')
grid()
```

We'll let $\dot{\xi}_{eq}(T)$ denote equilibrium value of the extent as a function of temperature. Those values are defined as roots to the equation

$$ \ln K_a(\dot{\xi}_{eq}(T)) - \ln K_r(T) = 0$$This implemented as a Python lambda function where a root=finding algorithm is used to solve the equilibrium condition as a function of temperature.

In [15]:

```
from scipy.optimize import brentq as fzero
xEquil = lambda T: fzero(lambda x: lnKa(x) - lnKr(T), 0, xMax)
```

Plotting

In [16]:

```
plot(T,map(xEquil,T))
xlabel('Temperature [K]')
ylabel('Extent of Reaction [gmol/min]')
title('Equilibrium Extent of Reaction at {:.1f} atm'.format(P))
grid()
```

In [17]:

```
yEquil = dict()
for c in C:
yEquil[c] = lambda T,c=c: y[c](xEquil(T))
for c in C:
plot(T,map(yEquil[c],T))
legend(C,'best')
xlabel('Temperature [K]')
ylabel('Mole Fraction [mol/mol]')
title('Reactor Outlet Gas Composition at {:.1f} atm'.format(P))
grid()
```

$$f_{conv} = \frac{\dot{n}_{in,CH_4} - \dot{n}_{out,CH_4}}{\dot{n}_{in,CH_4}}$$

In [18]:

```
fconv = lambda T: (nIn['CH4'] - nOut['CH4'](xEquil(T)))/nIn['CH4']
plot(T,map(fconv,T))
xlabel('Temperature [K]')
ylabel('Fractional Methane Conversion')
title('Equilibrium Conversion at {:.1f} atm'.format(P))
grid()
```

- Repeat the calculations for different operating pressures. How does equilibrium extent of reaction depend on operating pressure? Why is this the case?
- Repeat the calculations for different feed ratios for water and methane. How does extent of reaction depend on the feed ratio of water?
- Find component heat capacity data, then write a function to compute the heat requirement for the reactor. If the reactor is operated adiabatically, what is the temperature difference between the inlet and outlet streams?

In [37]:

```
```