## Method of continuity for nonlinear equation solving

| categories: nonlinear algebra | tags: | View Comments

Matlab post Adapted from Perry's Chemical Engineers Handbook, 6th edition 2-63.

We seek the solution to the following nonlinear equations:

$$2 + x + y - x^2 + 8 x y + y^3 = 0$$

$$1 + 2x - 3y + x^2 + xy - y e^x = 0$$

In principle this is easy, we simply need some initial guesses and a nonlinear solver. The challenge here is what would you guess? There could be many solutions. The equations are implicit, so it is not easy to graph them, but let us give it a shot, starting on the x range -5 to 5. The idea is set a value for x, and then solve for y in each equation.

import numpy as np
from scipy.optimize import fsolve

import matplotlib.pyplot as plt

def f(x, y):
return 2 + x + y - x**2 + 8*x*y + y**3;

def g(x, y):
return 1 + 2*x - 3*y + x**2 + x*y - y*np.exp(x)

x = np.linspace(-5, 5, 500)

@np.vectorize
def fy(x):
x0 = 0.0
def tmp(y):
return f(x, y)
y1, = fsolve(tmp, x0)
return y1

@np.vectorize
def gy(x):
x0 = 0.0
def tmp(y):
return g(x, y)
y1, = fsolve(tmp, x0)
return y1

plt.plot(x, fy(x), x, gy(x))
plt.xlabel('x')
plt.ylabel('y')
plt.legend(['fy', 'gy'])
plt.savefig('images/continuation-1.png')

>>> >>> >>> >>> ... ... >>> ... ... >>> >>> >>> ... ... ... ... ... ... ... >>> ... ... ... ... ... ... ... >>> >>> /opt/kitchingroup/enthought/epd-7.3-2-rh5-x86_64/lib/python2.7/site-packages/scipy/optimize/minpack.py:152: RuntimeWarning: The iteration is not making good progress, as measured by the
improvement from the last ten iterations.
warnings.warn(msg, RuntimeWarning)
/opt/kitchingroup/enthought/epd-7.3-2-rh5-x86_64/lib/python2.7/site-packages/scipy/optimize/minpack.py:152: RuntimeWarning: The iteration is not making good progress, as measured by the
improvement from the last five Jacobian evaluations.
warnings.warn(msg, RuntimeWarning)
[<matplotlib.lines.Line2D object at 0x1a0c4990>, <matplotlib.lines.Line2D object at 0x1a0c4a90>]
<matplotlib.text.Text object at 0x19d5e390>
<matplotlib.text.Text object at 0x19d61d90>
<matplotlib.legend.Legend object at 0x189df850>


You can see there is a solution near x = -1, y = 0, because both functions equal zero there. We can even use that guess with fsolve. It is disappointly easy! But, keep in mind that in 3 or more dimensions, you cannot perform this visualization, and another method could be required.

def func(X):
x,y = X
return [f(x, y), g(x, y)]

print fsolve(func, [-2, -2])

... ... >>> [ -1.00000000e+00   1.28730858e-15]


We explore a method that bypasses this problem today. The principle is to introduce a new variable, $$\lambda$$, which will vary from 0 to 1. at $$\lambda=0$$ we will have a simpler equation, preferrably a linear one, which can be easily solved, or which can be analytically solved. At $$\lambda=1$$, we have the original equations. Then, we create a system of differential equations that start at the easy solution, and integrate from $$\lambda=0$$ to $$\lambda=1$$, to recover the final solution.

We rewrite the equations as:

$$f(x,y) = (2 + x + y) + \lambda(- x^2 + 8 x y + y^3) = 0$$

$$g(x,y) = (1 + 2x - 3y) + \lambda(x^2 + xy - y e^x) = 0$$

Now, at $$\lambda=0$$ we have the simple linear equations:

$$x + y = -2$$

$$2x - 3y = -1$$

These equations are trivial to solve:

x0 = np.linalg.solve([[1., 1.], [2., -3.]],[ -2, -1])
print x0

[-1.4 -0.6]


We form the system of ODEs by differentiating the new equations with respect to $$\lambda$$. Why do we do that? The solution, (x,y) will be a function of $$\lambda$$. From calculus, you can show that:

$$\frac{\partial f}{\partial x}\frac{\partial x}{\partial \lambda}+\frac{\partial f}{\partial y}\frac{\partial y}{\partial \lambda}=-\frac{\partial f}{\partial \lambda}$$

$$\frac{\partial g}{\partial x}\frac{\partial x}{\partial \lambda}+\frac{\partial g}{\partial y}\frac{\partial y}{\partial \lambda}=-\frac{\partial g}{\partial \lambda}$$

Now, solve this for $$\frac{\partial x}{\partial \lambda}$$ and $$\frac{\partial y}{\partial \lambda}$$. You can use Cramer's rule to solve for these to yield:

\begin{eqnarray} \ \frac{\partial x}{\partial \lambda} &=& \frac{\partial f/\partial y \partial g/\partial \lambda - \partial f/\partial \lambda \partial g/\partial y}{\partial f/\partial x \partial g/\partial y - \partial f/\partial y \partial g/\partial x } \\\\ \frac{\partial y}{\partial \lambda} &=& \frac{\partial f/\partial \lambda \partial g/\partial x - \partial f/\partial x \partial g/\partial \lambda}{\partial f/\partial x \partial g/\partial y - \partial f/\partial y \partial g/\partial x } \end{eqnarray} For this set of equations: \begin{eqnarray} \ \partial f/\partial x &=& 1 - 2\lambda x + 8\lambda y \\\\ \partial f/\partial y &=& 1 + 8 \lambda x + 3 \lambda y^2 \\\\ \partial g/\partial x &=& 2 + 2 \lambda x + \lambda y - \lambda y e^x\\\\ \partial g/\partial y &=& -3 + \lambda x - \lambda e^x \end{eqnarray}

Now, we simply set up those two differential equations on $$\frac{\partial x}{\partial \lambda}$$ and $$\frac{\partial y}{\partial \lambda}$$, with the initial conditions at $$\lambda = 0$$ which is the solution of the simpler linear equations, and integrate to $$\lambda = 1$$, which is the final solution of the original equations!

def ode(X, LAMBDA):
x,y = X
pfpx = 1.0 - 2.0 * LAMBDA * x + 8 * LAMBDA * y
pfpy = 1.0 + 8.0 * LAMBDA * x + 3.0 * LAMBDA * y**2
pfpLAMBDA = -x**2 + 8.0 * x * y + y**3;
pgpx = 2. + 2. * LAMBDA * x + LAMBDA * y - LAMBDA * y * np.exp(x)
pgpy = -3. + LAMBDA * x - LAMBDA * np.exp(x)
pgpLAMBDA = x**2 + x * y - y * np.exp(x);
dxdLAMBDA = (pfpy * pgpLAMBDA - pfpLAMBDA * pgpy) / (pfpx * pgpy - pfpy * pgpx)
dydLAMBDA = (pfpLAMBDA * pgpx - pfpx * pgpLAMBDA) / (pfpx * pgpy - pfpy * pgpx)
dXdLAMBDA = [dxdLAMBDA, dydLAMBDA]
return dXdLAMBDA

from scipy.integrate import odeint

lambda_span = np.linspace(0, 1, 100)

X = odeint(ode, x0, lambda_span)

xsol, ysol = X[-1]
print 'The solution is at x={0:1.3f}, y={1:1.3f}'.format(xsol, ysol)
print f(xsol, ysol), g(xsol, ysol)

... ... ... ... ... ... ... ... ... >>> >>> >>> >>> >>> >>> >>> The solution is at x=-1.000, y=0.000
-1.27746598808e-06 -1.15873819107e-06


You can see the solution is somewhat approximate; the true solution is x = -1, y = 0. The approximation could be improved by lowering the tolerance on the ODE solver. The functions evaluate to a small number, close to zero. You have to apply some judgment to determine if that is sufficiently accurate. For instance if the units on that answer are kilometers, but you need an answer accurate to a millimeter, this may not be accurate enough.

This is a fair amount of work to get a solution! The idea is to solve a simple problem, and then gradually turn on the hard part by the lambda parameter. What happens if there are multiple solutions? The answer you finally get will depend on your $$\lambda=0$$ starting point, so it is possible to miss solutions this way. For problems with lots of variables, this would be a good approach if you can identify the easy problem.

org-mode source

## Calculating a bubble point pressure of a mixture

| categories: nonlinear algebra | tags: thermodynamics | View Comments

We previously learned to read a datafile containing lots of Antoine coefficients into a database, and use the coefficients to estimate vapor pressure of a single compound. Here we use those coefficents to compute a bubble point pressure of a mixture.

The bubble point is the temperature at which the sum of the component vapor pressures is equal to the the total pressure. This is where a bubble of vapor will first start forming, and the mixture starts to boil.

Consider an equimolar mixture of benzene, toluene, chloroform, acetone and methanol. Compute the bubble point at 760 mmHg, and the gas phase composition. The gas phase composition is given by: $$y_i = x_i*P_i/P_T$$.

import numpy as np
from scipy.optimize import fsolve

dtype=[('id', np.int),
('formula', 'S8'),
('name', 'S28'),
('A', np.float),
('B', np.float),
('C', np.float),
('Tmin', np.float),
('Tmax', np.float),
('??', 'S4'),
('?', 'S4')],
skiprows=7)

compounds = ['benzene', 'toluene', 'chloroform', 'acetone', 'methanol']

# extract the data we want
A = np.array([data[data['name'] == x]['A'][0] for x in compounds])
B = np.array([data[data['name'] == x]['B'][0] for x in compounds])
C = np.array([data[data['name'] == x]['C'][0] for x in compounds])
Tmin = np.array([data[data['name'] == x]['Tmin'][0] for x in compounds])
Tmax = np.array([data[data['name'] == x]['Tmax'][0] for x in compounds])

# we have an equimolar mixture
x = np.array([0.2, 0.2, 0.2, 0.2, 0.2])

# Given a T, we can compute the pressure of each species like this:

T = 67 # degC
P = 10**(A - B / (T + C))
print P
print np.dot(x, P)  # total mole-fraction weighted pressure

Tguess = 67
Ptotal = 760

def func(T):
P = 10**(A - B / (T + C))
return Ptotal - np.dot(x, P)

Tbubble, = fsolve(func, Tguess)

print 'The bubble point is {0:1.2f} degC'.format(Tbubble)

# double check answer is in a valid T range
if np.any(Tbubble < Tmin) or np.any(Tbubble > Tmax):
print 'T_bubble is out of range!'

# print gas phase composition
y = x * 10**(A - B / (Tbubble + C))/Ptotal

for cmpd, yi in zip(compounds, y):
print 'y_{0:<10s} = {1:1.3f}'.format(cmpd, yi)

[  498.4320267    182.16010994   898.31061294  1081.48181768   837.88860027]
699.654633507
The bubble point is 69.46 degC
y_benzene    = 0.142
y_toluene    = 0.053
y_chloroform = 0.255
y_acetone    = 0.308
y_methanol   = 0.242


org-mode source

## Solving CSTR design equations

| categories: nonlinear algebra | tags: reaction engineering | View Comments

Given a continuously stirred tank reactor with a volume of 66,000 dm^3 where the reaction $$A \rightarrow B$$ occurs, at a rate of $$-r_A = k C_A^2$$ ($$k=3$$ L/mol/h), with an entering molar flow of F_{A0} = 5 mol/h and a volumetric flowrate of 10 L/h, what is the exit concentration of A?

From a mole balance we know that at steady state $$0 = F_{A0} - F_A + V r_A$$. That equation simply states the sum of the molar flow of A in in minus the molar flow of A out plus the molar rate A is generated is equal to zero at steady state. This is directly the equation we need to solve. We need the following relationship:

1. $$F_A = v0 C_A$$
from scipy.optimize import fsolve

Fa0 = 5.0
v0 = 10.

V = 66000.0  # reactor volume L^3
k = 3.0      # rate constant L/mol/h

def func(Ca):
"Mole balance for a CSTR. Solve this equation for func(Ca)=0"
Fa = v0 * Ca     # exit molar flow of A
ra = -k * Ca**2  # rate of reaction of A L/mol/h
return Fa0 - Fa + V * ra

# CA guess that that 90 % is reacted away
CA_guess = 0.1 * Fa0 / v0
CA_sol, = fsolve(func, CA_guess)

print 'The exit concentration is {0} mol/L'.format(CA_sol)

The exit concentration is 0.005 mol/L


It is a little confusing why it is necessary to put a comma after the CA_sol in the fsolve command. If you do not put it there, you get brackets around the answer.

org-mode source

## Computing a pipe diameter

| categories: nonlinear algebra | tags: fluids | View Comments

Matlab post A heat exchanger must handle 2.5 L/s of water through a smooth pipe with length of 100 m. The pressure drop cannot exceed 103 kPa at 25 degC. Compute the minimum pipe diameter required for this application.

Adapted from problem 8.8 in Problem solving in chemical and Biochemical Engineering with Polymath, Excel, and Matlab. page 303.

We need to estimate the Fanning friction factor for these conditions so we can estimate the frictional losses that result in a pressure drop for a uniform, circular pipe. The frictional forces are given by $$F_f = 2f_F \frac{\Delta L v^2}{D}$$, and the corresponding pressure drop is given by $$\Delta P = \rho F_f$$. In these equations, $$\rho$$ is the fluid density, $$v$$ is the fluid velocity, $$D$$ is the pipe diameter, and $$f_F$$ is the Fanning friction factor. The average fluid velocity is given by $$v = \frac{q}{\pi D^2/4}$$.

For laminar flow, we estimate $$f_F = 16/Re$$, which is a linear equation, and for turbulent flow ($$Re > 2100$$) we have the implicit equation $$\frac{1}{\sqrt{f_F}}=4.0 \log(Re \sqrt{f_F})-0.4$$. Of course, we define $$Re = \frac{D v\rho}{\mu}$$ where $$\mu$$ is the viscosity of the fluid.

It is known that $$\rho(T) = 46.048 + 9.418 T -0.0329 T^2 +4.882\times10^{-5}-2.895\times10^{-8}T^4$$ and $$\mu = \exp\left({-10.547 + \frac{541.69}{T-144.53}}\right)$$ where $$\rho$$ is in kg/m^3 and $$\mu$$ is in kg/(m*s).

The aim is to find $$D$$ that solves: $$\Delta p = \rho 2 f_F \frac{\Delta L v^2}{D}$$. This is a nonlinear equation in $$D$$, since D affects the fluid velocity, the Re, and the Fanning friction factor. Here is the solution

import numpy as np
from scipy.optimize import fsolve
import matplotlib.pyplot as plt

T = 25 + 273.15
Q = 2.5e-3       # m^3/s
deltaP = 103000  # Pa
deltaL = 100     # m

#Note these correlations expect dimensionless T, where the magnitude
# of T is in K

def rho(T):
return 46.048 + 9.418 * T -0.0329 * T**2 +4.882e-5 * T**3 - 2.895e-8 * T**4

def mu(T):
return np.exp(-10.547 + 541.69 / (T - 144.53))

def fanning_friction_factor_(Re):
if Re < 2100:
raise Exception('Flow is probably not turbulent, so this correlation is not appropriate.')
# solve the Nikuradse correlation to get the friction factor
def fz(f): return 1.0/np.sqrt(f) - (4.0*np.log10(Re*np.sqrt(f))-0.4)
sol, = fsolve(fz, 0.01)
return sol

fanning_friction_factor = np.vectorize(fanning_friction_factor_)

Re = np.linspace(2200, 9000)
f = fanning_friction_factor(Re)

plt.plot(Re, f)
plt.xlabel('Re')
plt.ylabel('fanning friction factor')
# You can see why we use 0.01 as an initial guess for solving for the
# Fanning friction factor; it falls in the middle of ranges possible
# for these Re numbers.
plt.savefig('images/pipe-diameter-1.png')

def objective(D):
v = Q / (np.pi * D**2 / 4)
Re = D * v * rho(T) / mu(T)

fF = fanning_friction_factor(Re)

return deltaP - 2 * fF * rho(T) * deltaL * v**2 / D

D, = fsolve(objective, 0.04)

print('The minimum pipe diameter is {0} m\n'.format(D))

The minimum pipe diameter is 0.0389653369531 m


Any pipe diameter smaller than that value will result in a larger pressure drop at the same volumetric flow rate, or a smaller volumetric flowrate at the same pressure drop. Either way, it will not meet the design specification.

org-mode source

## Water gas shift equilibria via the NIST Webbook

| categories: nonlinear algebra | tags: | View Comments

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:

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 H_2 at 1000K.

import numpy as np

T = np.linspace(500,1000) # degrees K
t = T/1000;


## 1 hydrogen

# 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*np.log(t) + B*t + C*t**2/2 + D*t**3/3 - E/(2*t**2) + G);


## 2 H_{2}O

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*np.log(t) + B*t + C*t**2/2 + D*t**3/3 - E/(2*t**2) + G);


## 3 CO

# 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*np.log(t) + B*t + C*t**2/2 + D*t**3/3 - E/(2*t**2) + G);


## 4 CO_{2}

# 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*np.log(t) + B*t + C*t**2/2 + D*t**3/3 - E/(2*t**2) + G);


## 5 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;

print('deltaH = {0:1.2f}'.format(Hrxn_29815))
print('deltaG = {0:1.2f}'.format(Grxn_29815))

>>> >>> >>> deltaH = -41.15
deltaG = -28.62


## 6 Non-standard state $$\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

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.

Hrxn = Hrxn_29815 + dH_CO2 + dH_H2 - dH_CO - dH_H2O
Grxn = Hrxn - T*(S_CO2 + S_H2 - S_CO - S_H2O)/1000


## 7 Plot how the $$\Delta G$$ varies with temperature

import matplotlib.pyplot as plt
plt.figure(); plt.clf()
plt.plot(T,Grxn, label='$\Delta G_{rxn}$')
plt.plot(T,Hrxn, label='$\Delta H_{rxn}$')
plt.xlabel('Temperature (K)')
plt.ylabel('(kJ/mol)')
plt.legend( loc='best')
plt.savefig('images/wgs-nist-1.png')

<matplotlib.figure.Figure object at 0x04199CF0>
[<matplotlib.lines.Line2D object at 0x0429BF30>]
[<matplotlib.lines.Line2D object at 0x0427DFB0>]
<matplotlib.text.Text object at 0x041B79F0>
<matplotlib.text.Text object at 0x040CEF70>
<matplotlib.legend.Legend object at 0x043CB5F0>


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.

## 8 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 = np.exp(-Grxn/R/T);

plt.figure()
plt.plot(T,K)
plt.xlim([500, 1000])
plt.xlabel('Temperature (K)')
plt.ylabel('Equilibrium constant')
plt.savefig('images/wgs-nist-2.png')

>>> >>> <matplotlib.figure.Figure object at 0x044DBE90>
[<matplotlib.lines.Line2D object at 0x045A53F0>]
(500, 1000)
<matplotlib.text.Text object at 0x04577470>
<matplotlib.text.Text object at 0x0457F410>


## 9 Equilibrium yield of WGS

Now let us suppose we have a reactor with a feed of H_2O and CO at 10atm at 1000K. What is the equilibrium yield of H_2? 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.

from scipy.interpolate import interp1d
from scipy.optimize import fsolve

#
# 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 simply computed a vector of enthalpies, entropies, etc... Here we interpolate
K_func = interp1d(T,K);
K_Temperature = K_func(1000)

# 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.

def f(X):
return K_Temperature - X**2/(1-X)**2;

x0 = 0.5
Xeq, = fsolve(f, x0)

print('The equilibrium conversion for these feed conditions is: {0:1.2f}'.format(Xeq))

>>> >>> ... ... ... ... ... >>> >>> >>> >>> >>> ... ... ... >>> >>> >>> >>> ... ... ... ... >>> ... ... >>> >>> >>> 0.54504291144
The equilibrium conversion for these feed conditions is: 0.55


## 10 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 H_2 and CO_2 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 H_2 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

print P_CO,P_H2O, P_H2, P_CO2

>>> >>> >>> >>> 2.2747854428 2.2747854428 2.7252145572 2.7252145572


## 11 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!

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

1.43522674762
1.43522674762


They are the same.

## 12 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.

I noticed a different behavior in interpolation between scipy.interpolate.interp1d and Matlab's interp1. The scipy function returns an interpolating function, whereas the Matlab function directly interpolates new values, and returns the actual interpolated data.