This notebook shows how to solve chemical kintics problems for a continuously stirred tank reactor using ChemPy.
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
rsys = ReactionSystem.from_string("A -> B; 'k'", substance_factory=lambda k: Substance(k))
rsys
odesys, extra = get_odesys(rsys, include_params=False)
odesys.names, odesys.param_names
We can change the expressions of the ODE system manually to account for source and sink terms from the flux:
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)))) + '$')
cstr.param_names
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')
res.plot()
We can derive an analytic solution to the ODE system:
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
exprs0 = [expr.subs(t, 0) for expr in exprs]
exprs0
sol = cstr.be.solve([expr - c0 for expr, c0 in zip(exprs0, (IA, IB))], (c1, c2))
sol
exprs2 = [expr.subs(sol) for expr in exprs]
exprs2
IA
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 ('))
exprs2_0 = [expr.subs(t, 0).simplify() for expr in exprs2]
exprs2_0
_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'])
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]
yref = get_ref(res)
yref.shape
Plotting the error (by comparison to the analytic solution):
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)
ChemPy has support for generating the CSTR expressions directly in ReactionSystem.rates
& get_odesys
. This simplifies the code the user needs to write considerably:
cstr2, extra2 = get_odesys(rsys, include_params=False, cstr=True)
cstr2.exprs
Note how we only needed to pass cstr=True
to get_odesys
to get the right expressions.
cstr2.param_names
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)
The analytic solution of a unary reaction in a CSTR is also available in ChemPy:
from chempy.kinetics.integrated import unary_irrev_cstr
help(unary_irrev_cstr)
The symbols of the feedratio and feed concentrations are available in the second output of get_odesys
(the extra dictionary):
fr, fc = extra2['cstr_fr_fc']
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]
res2.plot(y=res2.yout - get_ref2(res2))
assert np.allclose(res2.yout, get_ref2(res2))