In [ ]:
import sympy as sp
sp.init_printing()
Analyitc solution for CSTR of 2 A -> B
In [ ]:
symbs = t, f, fcA, fcB, IA, IB, k, c1, c2 = sp.symbols('t f phi_A phi_B I_A I_B k c1 c2', real=True, positive=True)
symbs
In [ ]:
def analytic(t, f, fcA, fcB, IA, IB, k, c1, c2):
u = sp.sqrt(f*(f + 4*fcA*k))
v = u*sp.tanh(u/2*(t - c1))
return [
(-f + v)/(2*k),
sp.exp(-f*t)*c2 + (f + 2*(fcA + fcB)*k - v)/(2*k)
]
In [ ]:
exprs = analytic(*symbs)
exprs
In [ ]:
exprs0 = [expr.subs(t, 0) for expr in exprs]
exprs0
In [ ]:
sol = sp.solve([expr - c0 for expr, c0 in zip(exprs0, [IA, IB])], [c1, c2])
sol
In [ ]:
exprs2 = [expr.subs(dict(zip([c1, c2], sol[1]))) for expr in exprs]
exprs2
In [ ]:
[expr.subs(t, 0).simplify() for expr in exprs2]
In [ ]:
%matplotlib inline
import numpy as np
from chempy import ReactionSystem
from chempy.kinetics.ode import get_odesys
rsys = ReactionSystem.from_string("OH + OH -> H2O2; 'k'")
cstr, extra = get_odesys(rsys, include_params=False, cstr=True)
fr, fc = extra['cstr_fr_fc']
print(cstr.names, cstr.param_names)
cstr.exprs
In [ ]:
args = tout, c0, params = np.linspace(0, .17), {'OH': 2, 'H2O2': 3}, {'k': 5, fc['OH']: 42, fc['H2O2']: 11, fr: 13}
res = cstr.integrate(*args)
res.plot()
In [ ]:
def analytic_alt(t, f, fcA, fcB, IA, IB, k):
u = sp.sqrt(f*(f + 4*fcA*k))
q = sp.atanh(-u*(sp.sqrt(f) + 2*k*IA)/(f*(f + 4*fcA*k)))
v = u*sp.tanh(u/2*t - q)
w = sp.exp(-f*t)/2/k
y = 2*k*(fcA + fcB)
return [
(-f + v)/(2*k),
w*(sp.exp(f*t)*(f + y - v) - y + 2*k*(IA + IB))
]
In [ ]:
def analytic_alt2(t, f, fcA, fcB, IA, IB, k, n):
one_point_five = sp.S(3)/2
a0, b0, x, y = fcA, fcB, IA, IB
Sqrt, Tanh, ArcTanh, E = sp.sqrt, sp.tanh, sp.atanh, sp.E
return [
(-f + Sqrt(f)*Sqrt(f + 8*a0*k)*
Tanh((Sqrt(f)*Sqrt(f + 8*a0*k)*t - 2*ArcTanh((-(f**one_point_five*Sqrt(f + 8*a0*k)) - 4*Sqrt(f)*k*Sqrt(f + 8*a0*k)*x)/(f**2 + 8*a0*f*k)))/2)
)/(4*k),
(-8*b0*k + 8*b0*E**(f*t)*k + E**(f*t)*f*n - 4*a0*k*n + 4*a0*E**(f*t)*k*n + 4*k*n*x + 8*k*y -
E**(f*t)*Sqrt(f)*Sqrt(f + 8*a0*k)*n*Tanh((Sqrt(f)*Sqrt(f + 8*a0*k)*
(t - (2*ArcTanh((-(f**one_point_five*Sqrt(f + 8*a0*k)) - 4*Sqrt(f)*k*Sqrt(f + 8*a0*k)*x)/(f**2 + 8*a0*f*k)))/(Sqrt(f)*Sqrt(f + 8*a0*k))))
/2))/(8*E**(f*t)*k)
]
# return [
# (-f + Sqrt(f)*Sqrt(f + 8*a0*k)*
# Tanh((Sqrt(f)*Sqrt(f + 8*a0*k)*t - 2*ArcTanh((-(f**one_point_five*Sqrt(f + 8*a0*k)) - 4*Sqrt(f)*k*Sqrt(f + 8*a0*k)*x)/(f**2 + 8*a0*f*k)))/2)
# )/(4*k),
# (E**(f*t)*f - 4*a0*k - 8*b0*k + 4*a0*E**(f*t)*k + 8*b0*E**(f*t)*k + 4*k*x + 8*k*y -
# E**(f*t)*Sqrt(f)*Sqrt(f + 8*a0*k)*Tanh((Sqrt(f)*Sqrt(f + 8*a0*k)*
# (t - (2*ArcTanh((-(f**one_point_five*Sqrt(f + 8*a0*k)) - 4*Sqrt(f)*k*Sqrt(f + 8*a0*k)*x)/(f**2 + 8*a0*f*k)))/(Sqrt(f)*Sqrt(f + 8*a0*k))))
# /2))/(8*E**(f*t)*k)
# ]
In [ ]:
n = sp.Symbol('n')
exprs_alt = analytic_alt2(*symbs[:-2], n)
exprs_alt
In [ ]:
cses, expr_cse = sp.cse([expr.subs({fcA: sp.Symbol('fr'), fcB: sp.Symbol('fp'), f: sp.Symbol('fv'),
IA: sp.Symbol('r'), IB: sp.Symbol('p')}) for expr in exprs_alt])
In [ ]:
print(
'\n'.join(['%s = %s' % (lhs, rhs) for lhs, rhs in cses] + ['return (\n %s\n)' % str(expr_cse)[1:-1]]).replace(
'3/2', 'three/2').replace(
'exp', 'be.exp').replace(
'sqrt', 'be.sqrt').replace(
'atanh', 'ATANH').replace(
'tanh', 'be.tanh\n ').replace(
'ATANH', 'atanh')
)
In [ ]:
[expr.subs(t, 0).simplify() for expr in exprs_alt]
In [ ]:
[expr.diff(t).subs(t, 0).simplify() for expr in exprs_alt]
In [ ]:
print(list(rsys.substances))
rsys
In [ ]:
print(cstr.names, cstr.param_names)
cstr.exprs
In [ ]:
symbs[:-2]
In [ ]:
_cb = sp.lambdify(symbs[:-2] + (n,), exprs_alt)
def calc_analytic(xout, y0, p):
return _cb(xout, p[fr], p[fc['OH']], p[fc['H2O2']], y0['OH'], y0['H2O2'], p['k'], 1)
In [ ]:
def get_analytic(result):
ref = calc_analytic(
result.xout,
{k: res.named_dep(k)[0] for k in result.odesys.names},
{k: res.named_param(k) for k in result.odesys.param_names})
return np.array([ref[{'OH': 0, 'H2O2': 1}[k]] for k in result.odesys.names]).T
In [ ]:
yref = get_analytic(res)
print(yref.shape)
In [ ]:
res.plot()
res.plot(y=yref)
In [ ]:
res.plot(y=res.yout - yref)
In [ ]:
from chempy.kinetics.integrated import binary_irrev_cstr
In [ ]:
def get_analytic2(result):
ref = binary_irrev_cstr(result.xout, result.named_param('k'), result.named_dep('OH')[0],
result.named_dep('H2O2')[0], result.named_param(fc['OH']), result.named_param(fc['H2O2']),
result.named_param(fr))
return np.array([ref[{'OH': 0, 'H2O2': 1}[k]] for k in result.odesys.names]).T
In [ ]:
res.plot(y=res.yout - get_analytic2(res))
In [ ]:
# Pseduo reversible reaction (batch reactor)
#
# A + B <-> C
# t a0 b0-x c0 + x
import sympy
kf, kb, a0, b0, c0, x = sympy.symbols('k_{\\rm\\ f} k_{\\rm\\ b} a0 b0 c0 x', real=True, nonnegative=True)
rf = kf*a0*(b0-x)
rb = kb*(c0+x)
dxdt = (rf - rb).expand().factor(x)
dxdt
In [ ]:
integral = sympy.integrate(1/dxdt, (x, 0, x))
integral
In [ ]:
sol, = sympy.solve(integral - t, x)
sol
In [ ]:
print(sympy.ccode(sol.subs({kf: 'kf', kb: 'kb', a0: 'major', b0: 'minor', c0: 'prod'})))