[Python] How to apply the method of lines (MOL) with scipy odeint function with Danckwerts...

I am trying to solve a convection-diffusion-reaction problem for a tubular reactor. The reaction governing the process is as follows:



The method I chose to solve it is the line method, which consists of solving a system of ordinary differential equations in time, with discretization in space. I found a topic here solving-system-of-des-using-method-of-lines-spatially-discretised that deals with a similar problem, but the biggest difficulty I am encountering is related to my boundary conditions (Danckwerts):



From the first boundary condition, the value of C1 can be obtained with a second-order progressive approximation:



From the second, using a second-order regressive approximation, the value of C_N+1 can be obtained:



I tried the following code, but from what I read in the scipy documentation scipy.integrate.odeint for the odeint function, it is apparently not correct to modify the y array within the function, but only its derivatives.

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm
from matplotlib.ticker import LinearLocator
import time

# Parameters
L = 1.0 # meter (reator length)
v_x = 3.3e-4 # m/s (average fluid velocity)
D = 0.0000033 # m^2/s (diffusion coefficient)
Pe = L*v_x/D # Peclet number

kd = 0.000055 # 1/s (dissociation constant)
factor = 0.6 # efficient factor
Da = L*kd*factor/v_x # Damköhler number

kt = 6e7
Da_mu0 = L*kt/v_x

# Grid and steps
m = 101
n = 100001
x = np.linspace(0, 1, m) # spatial
t = np.linspace(0, 1, n) # temporal
h = x[1] - x[0] # spatial step
k = t[1] - t[0] # temporal step

from scipy.integrate import odeint

init_I = np.ones([m])
init_mu0 = np.zeros([m])

def odefunc(y,t):

I, mu0 = np.split(y,2)

dIdt = np.zeros(x.shape)
dmu0dt = np.zeros(x.shape)

I[0] = (2*Pe/m + 4*I[1] - I[2])/(2*Pe/m + 3)
I[-1] = (4*I[-2] - I[-3])/3

mu0[0] = (4*mu0[1] - mu0[2])/(3 + 2*Pe/m)
mu0[-1] = (4*mu0[-2] - mu0[-3])/3

for i in range(1, m-1):
dIdt = -(I[i+1] - I[i-1])/(2/m) + (1/Pe)*(I[i+1] - 2*I + I[i-1])/(1/m**2) - Da*(I[i+1] + I[i-1])/2
dmu0dt = -(mu0[i+1] - mu0[i-1])/(2/m) + (1/Pe)*(mu0[i+1] - 2*mu0 + mu0[i-1])/(1/m**2) + Da*(I[i+1] + I[i-1]) - (Da_mu0/4)*(mu0[i+1] +mu0[i-1])**2

return np.ravel([dIdt, dmu0dt])

init = np.ravel([init_I, init_mu0])
sol = odeint(odefunc, init, t, atol=1e-10, rtol=1e-10)
sol_I = sol[:,:m]
sol_mu0 = sol[:,m:]

fig, ax = plt.subplots(subplot_kw={"projection": "3d"})
M, N = np.meshgrid(x, t)

surf = ax.plot_surface(M, N, sol_I, cmap=cm.viridis, linewidth=0, antialiased=True)

ax.set_zlim(0, 1.0)
ax.zaxis.set_major_locator(LinearLocator(5))
ax.set_xlabel('x\'')
ax.set_ylabel('t\'')
ax.set_zlabel('c\'')

ax.zaxis.set_major_formatter('{x:.02f}')
ax.set_title(f'{Pe =:.02f} and {Da =:.02f}')

fig.colorbar(surf, shrink=0.5, aspect=5)

plt.show()

fig, ax = plt.subplots(subplot_kw={"projection": "3d"})

surf = ax.plot_surface(M, N, sol_mu0, cmap=cm.viridis, linewidth=0, antialiased=True)

ax.set_zlim(0, 1.0)
ax.zaxis.set_major_locator(LinearLocator(5))
ax.set_xlabel('x\'')
ax.set_ylabel('t\'')
ax.set_zlabel('c\'')

ax.zaxis.set_major_formatter('{x:.02f}')
ax.set_title(f'{Pe =:.02f} and {Da =:.02f}')

fig.colorbar(surf, shrink=0.5, aspect=5)

plt.show()

Continue reading...