Continuously stirred tank reactor (CSTR)

This notebook shows how to solve chemical kintics problems for a continuously stirred tank reactor using ChemPy.

In [1]:
from collections import defaultdict
import numpy as np
from IPython.display import Latex
import matplotlib.pyplot as plt
from pyodesys.symbolic import SymbolicSys
from chempy import Substance, ReactionSystem
from chempy.kinetics.ode import get_odesys
from chempy.units import SI_base_registry, default_units as u
from chempy.util.graph import rsys2graph
%matplotlib inline
In [2]:
rsys = ReactionSystem.from_string("A -> B; 'k'", substance_factory=lambda k: Substance(k))
rsys
Out[2]:
AB MassAction((Symbol((), ('k',)),))
In [3]:
odesys, extra = get_odesys(rsys, include_params=False)
odesys.names, odesys.param_names
Out[3]:
(('A', 'B'), ('k',))

We can change the expressions of the ODE system manually to account for source and sink terms from the flux:

In [4]:
t, c1, c2, IA, IB, f, k, fc_A, fc_B = map(odesys.be.Symbol, 't c1 c2 I_A I_B f k phi_A phi_B'.split())
newp = f, fc_A, fc_B
c_feed = {'A': fc_A, 'B': fc_B}
cstr = SymbolicSys(
    [
        (dep, expr - f*dep + f*c_feed[name]) for name, dep, expr
        in zip(odesys.names, odesys.dep, odesys.exprs)
    ],
    params=list(odesys.params) + list(newp),
    names=odesys.names,
    param_names=list(odesys.param_names) + [p.name for p in newp],
    dep_by_name=True,
    par_by_name=True,
)
Latex('$' + r'\\ '.join(map(lambda x: ':'.join(x), zip(map(lambda x: r'\frac{d%s}{dt}' % x, cstr.names),
                                                       map(cstr.be.latex, cstr.exprs)))) + '$')
Out[4]:
$\frac{dA}{dt}:f \phi_{A} - f y_{0} - p_{0} y_{0}\\ \frac{dB}{dt}:f \phi_{B} - f y_{1} + p_{0} y_{0}$
In [5]:
cstr.param_names
Out[5]:
('k', 'f', 'phi_A', 'phi_B')
In [6]:
init_c, pars = {'A': .15, 'B': .1}, {'k': 0.8, 'f': 0.3, 'phi_A': .7, 'phi_B': .1}
res = cstr.integrate(10, init_c, pars, integrator='cvode')
In [7]:
res.plot()
Out[7]:
<matplotlib.axes._subplots.AxesSubplot at 0x7ff0a35da5f8>

We can derive an analytic solution to the ODE system:

In [8]:
k = cstr.params[0]
e = cstr.be.exp
exprs = [
    fc_A*f/(f + k) + c1 * e(-t*(f + k)),
    (fc_A*k + fc_B*(f + k))/(f + k) - c1*e(-f*t)*(e(-t*k) - 1) + c2*e(-f*t)
]
cstr.be.init_printing()
exprs
Out[8]:
$\displaystyle \left[ c_{1} e^{- t \left(f + p_{0}\right)} + \frac{f \phi_{A}}{f + p_{0}}, \ - c_{1} \left(-1 + e^{- p_{0} t}\right) e^{- f t} + c_{2} e^{- f t} + \frac{p_{0} \phi_{A} + \phi_{B} \left(f + p_{0}\right)}{f + p_{0}}\right]$
In [9]:
exprs0 = [expr.subs(t, 0) for expr in exprs]
exprs0
Out[9]:
$\displaystyle \left[ c_{1} + \frac{f \phi_{A}}{f + p_{0}}, \ c_{2} + \frac{p_{0} \phi_{A} + \phi_{B} \left(f + p_{0}\right)}{f + p_{0}}\right]$
In [10]:
sol = cstr.be.solve([expr - c0 for expr, c0 in zip(exprs0, (IA, IB))], (c1, c2))
sol
Out[10]:
$\displaystyle \left\{ c_{1} : \frac{I_{A} f + I_{A} p_{0} - f \phi_{A}}{f + p_{0}}, \ c_{2} : \frac{I_{B} f + I_{B} p_{0} - f \phi_{B} - p_{0} \phi_{A} - p_{0} \phi_{B}}{f + p_{0}}\right\}$
In [11]:
exprs2 = [expr.subs(sol) for expr in exprs]
exprs2
Out[11]:
$\displaystyle \left[ \frac{f \phi_{A}}{f + p_{0}} + \frac{\left(I_{A} f + I_{A} p_{0} - f \phi_{A}\right) e^{- t \left(f + p_{0}\right)}}{f + p_{0}}, \ - \frac{\left(-1 + e^{- p_{0} t}\right) \left(I_{A} f + I_{A} p_{0} - f \phi_{A}\right) e^{- f t}}{f + p_{0}} + \frac{p_{0} \phi_{A} + \phi_{B} \left(f + p_{0}\right)}{f + p_{0}} + \frac{\left(I_{B} f + I_{B} p_{0} - f \phi_{B} - p_{0} \phi_{A} - p_{0} \phi_{B}\right) e^{- f t}}{f + p_{0}}\right]$
In [12]:
IA
Out[12]:
$\displaystyle I_{A}$
In [13]:
import sympy as sp
cses, expr_cse = sp.cse([expr.subs({fc_A: sp.Symbol('fr'), fc_B: sp.Symbol('fp'), f: sp.Symbol('fv'),
                                     IA: sp.Symbol('r'), IB: sp.Symbol('p')}) for expr in exprs2])
s = '\n'.join(['%s = %s' % (lhs, rhs) for lhs, rhs in cses] + [str(tuple(expr_cse))])
print(s.replace('exp', 'be.exp').replace('\n(', '\nreturn ('))
x0 = fv + p_0
x1 = 1/x0
x2 = fr*fv
x3 = x1*(fv*r + p_0*r - x2)
x4 = fr*p_0
x5 = be.exp(-fv*t)
return (x1*x2 + x3*be.exp(-t*x0), x1*x5*(-fp*fv - fp*p_0 + fv*p + p*p_0 - x4) + x1*(fp*x0 + x4) - x3*x5*(-1 + be.exp(-p_0*t)))
In [14]:
exprs2_0 = [expr.subs(t, 0).simplify() for expr in exprs2]
exprs2_0
Out[14]:
$\displaystyle \left[ I_{A}, \ I_{B}\right]$
In [15]:
_cb = cstr.be.lambdify([t, IA, IB, k, f, fc_A, fc_B], exprs2)
def analytic(x, c0, params):
    return _cb(x, c0['A'], c0['B'], params['k'], params['f'], params['phi_A'], params['phi_B'])
In [16]:
def get_ref(result, parameters=None):
    drctn = -1 if result.odesys.names[0] == 'B' else 1
    return np.array(analytic(
        result.xout,
        {k: result.named_dep(k)[0] for k in result.odesys.names},
        parameters or {k: result.named_param(k) for k in result.odesys.param_names})).T[:, ::drctn]
In [17]:
yref = get_ref(res)
yref.shape
Out[17]:
$\displaystyle \left( 90, \ 2\right)$

Plotting the error (by comparison to the analytic solution):

In [18]:
fig, axes = plt.subplots(1, 2, figsize=(14, 4))
res.plot(ax=axes[0])
res.plot(ax=axes[0], y=yref)
res.plot(ax=axes[1], y=res.yout - yref)
Out[18]:
<matplotlib.axes._subplots.AxesSubplot at 0x7ff0bfc8c748>

Automatically generating CSTR expressions using chempy

ChemPy has support for generating the CSTR expressions directly in ReactionSystem.rates & get_odesys. This simplifies the code the user needs to write considerably:

In [19]:
cstr2, extra2 = get_odesys(rsys, include_params=False, cstr=True)
cstr2.exprs
Out[19]:
$\displaystyle \left( p_{0} \left(p_{2} - y_{0}\right) - p_{3} y_{0}, \ p_{0} \left(p_{1} - y_{1}\right) + p_{3} y_{0}\right)$

Note how we only needed to pass cstr=True to get_odesys to get the right expressions.

In [20]:
cstr2.param_names
Out[20]:
('feedratio', 'fc_B', 'fc_A', 'k')
In [21]:
renamed_pars = {'k': pars['k'], 'fc_A': pars['phi_A'], 'fc_B': pars['phi_B'], 'feedratio': pars['f']}
res2 = cstr2.integrate(10, init_c, renamed_pars, integrator='cvode')
ref2 = get_ref(res2, pars)
res2.plot(y=res2.yout - ref2)
Out[21]:
<matplotlib.axes._subplots.AxesSubplot at 0x7ff0bfbeecc0>

The analytic solution of a unary reaction in a CSTR is also available in ChemPy:

In [22]:
from chempy.kinetics.integrated import unary_irrev_cstr
help(unary_irrev_cstr)
Help on function unary_irrev_cstr in module chempy.kinetics.integrated:

unary_irrev_cstr(t, k, r, p, fr, fp, fv, backend=None)
    Analytic solution for ``A -> B`` in a CSTR.
    
    Analytic solution for a first order process in a continuously
    stirred tank reactor (CSTR).
    
    Parameters
    ----------
    t : array_like
    k : float_like
        Rate constant
    r : float_like
        Initial concentration of reactant.
    p : float_like
        Initial concentration of product.
    fr : float_like
        Concentration of reactant in feed.
    fp : float_like
        Concentration of product in feed.
    fv : float_like
        Feed rate / tank volume ratio.
    backend : module or str
        Default is 'numpy', can also be e.g. ``sympy``.
    
    Returns
    -------
    length-2 tuple
        concentrations of reactant and product

The symbols of the feedratio and feed concentrations are available in the second output of get_odesys (the extra dictionary):

In [23]:
fr, fc = extra2['cstr_fr_fc']
In [24]:
def get_ref2(result):
    drctn = -1 if result.odesys.names[0] == 'B' else 1
    return np.array(unary_irrev_cstr(
        result.xout,
        result.named_param('k'),
        result.named_dep('A')[0],
        result.named_dep('B')[0], 
        result.named_param(fc['A']),
        result.named_param(fc['B']),
        result.named_param(fr))).T[:, ::drctn]
In [25]:
res2.plot(y=res2.yout - get_ref2(res2))
Out[25]:
<matplotlib.axes._subplots.AxesSubplot at 0x7ff0bfb4ea90>
In [26]:
assert np.allclose(res2.yout, get_ref2(res2))
In [ ]: