The image below illustrates the kinetic scheme I am trying to model.
My first example focuses on the conversion of the wood-oil to non-volatiles and volatiles. The rate equations are as follows:
$$ r_w = \frac{dC_w}{dt} = -(K_{nv} + K_v)\,C_w \\ r_{nv} = \frac{dC_{nv}}{dt} = K_{nv}\,C_w \\ r_{v} = \frac{dC_{v}}{dt} = K_{v}\,C_w \\ $$
To solve for the concentrations, I use the odeint
solver in SciPy as demonstrated in the Python code below.
dt = 0.0001 # time step, delta t
tmax = 25 # max time, s
t = np.linspace(0, tmax, num=tmax/dt) # time vector
nt = len(t) # total number of time steps
def rates(c, t):
"""
w = wood-oil as conc[0]
nv = non-volatiles as conc[1]
v = volatiles as conc[2]
"""
Knv = 0.3
Kv = 0.8
rw = -(Knv + Kv) * c[0]
rnv = Knv * c[0]
rv = Kv * c[0]
return [rw, rnv, rv]
cc = sp.odeint(rates, [1, 0, 0], t)
The above code displays the following plot.
In my next approach I attempt to account for the reactions that consume the non-volatiles and volatiles; the rate equations are displayed below.
$$ r_w = \frac{dC_w}{dt} = -(K_{nv} + K_v)\,C_w \\ r_{nv} = \frac{dC_{nv}}{dt} = K_{nv}\,C_w - K_{r1}\,{C_{nv}}^{2.5} - K_{cr}\,{C_{nv}}^{0.9} - K_{c1}\,{C_{nv}}^{1.1} \\ r_{v} = \frac{dC_{v}}{dt} = K_{v}\,C_w + K_{cr}\,{C_{nv}}^{0.9} - K_{d}\,{C_{v}}^{0.9} - K_{a}\,{C_{v}}^{1.4} - K_{g}\,{C_{v}}^{0.8} - K_{r2}\,{C_{v}}^{0.7} $$
And the code I use to solve for the concentrations is shown below.
def dCdt(c, t):
"""
w = wood-oil as conc[0]
nv = non-volatiles as conc[1]
v = volatiles as conc[2]
"""
Knv = 0.3
Kv = 0.8
Kr1 = 9.2e7; r1 = 2.5
Kcr = 4.1e-5; cr = 0.9
Kc1 = 3.7e5; c1 = 1.1
Kd = 8.0e-4; d = 0.9
Ka = 6.1e-6; a = 1.4
Kg = 1.8e-4; g = 0.8
Kr2 = 37.0e5; r2 = 0.7
rw = -(Knv + Kv)*c[0]
rnv = Knv*c[0] - Kr1*c[1]**r1 - Kcr*c[1]**cr - Kc1*c[1]**c1
rv = Kv*c[0] + Kcr*c[1]**cr - Kd*c[2]**d - Ka*c[2]**a - Kg*c[2]**g - Kr2*c[2]**r2
return [rw, rnv, rv]
cc2 = sp.odeint(dCdt, [1, 0, 0], t)
Unfortunately, this approach does not work and I get an error from Python about invalid value encountered in double_scalars
. It appears that the odeint
solver isn't capable of handling the different reaction orders.
So my next approach is to solve the system with the SciPy ode
solver. The code for this is shown below:
def dcdt(t, c):
"""
w = wood-oil as c[0]
nv = non-volatiles as c[1]
v = volatiles as c[2]
"""
Knv = 0.3
Kv = 0.8
Kr1 = 9.2e7; r1 = 2.5
Kcr = 4.1e-5; cr = 0.9
Kc1 = 3.7e5; c1 = 1.1
Kd = 8.0e-4; d = 0.9
Ka = 6.1e-6; a = 1.4
Kg = 1.8e-4; g = 0.8
Kr2 = 37.0e5; r2 = 0.7
rw = -(Knv + Kv)*c[0]
rnv = Knv*c[0] - Kr1*c[1]**r1 - Kcr*c[1]**cr - Kc1*c[1]**c1
rv = Kv*c[0] + Kcr*c[1]**cr - Kd*c[2]**d - Ka*c[2]**a - Kg*c[2]**g - Kr2*c[2]**r2
return [rw, rnv, rv]
# store concentrations
Coil = np.ones(nt) # bio-oil concentration
Cnvol = np.zeros(nt) # non-volatiles concentration
Cvol = np.zeros(nt) # volatiles concentration
# Setup the ode integrator where 'dopri5' is Runge-Kutta 4th order
r = sp.ode(dcdt).set_integrator('dopri5', nsteps=10000)
r.set_initial_value([1, 0, 0], 0)
# integrate the odes for each time step then store the results
k = 1
while r.successful() and r.t < tmax-dt:
r.integrate(r.t+dt)
Coil[k] = r.y[0]
Cnvol[k] = r.y[1]
Cvol[k] = r.y[2]
k += 1
Unfortunately, the ode
approach does not work and I receive a warning about the system being stiff. Does anyone have suggestions on how to solve this system of rate equations in Python when the reaction order is not one?
ode
solver you suggested but it will not converge to a solution. See my updated question for the new code. $\endgroup$ – wigging Jun 9 '16 at 16:42