In [36]:
import numpy as np
import pandas as pd
import math
import cmath
from scipy.optimize import root
from scipy.integrate import odeint
from __future__ import division
from scipy import *
from scipy.optimize import curve_fit
from pylab import *
import matplotlib.pyplot as plt
%matplotlib inline
import pylab as pp
from scipy import integrate, interpolate
from scipy import optimize
L'extraction supercritique est de plus en plus utilisée afin de retirer des matières organiques de différents liquides ou matrices solides. Cela est dû au fait que les fluides supercritiques ont des avantages non négligeables par rapport aux autres solvants, ils ont des caractèreistiques comprises entre celles des gaz et celles des solides. En changeant la température et la pression ils peuvent capter des composés différents, ils sont donc très efficaces. Le méchanisme de l'extraction supercritique est le suivant :
A - Le modèle de Reverchon :
Afin d'utiliser ce modèle, définissons les variables qui vont y être admises, ci-dessous la nomenclature du modèle :
Le modèle : Il est basé sur l'intégration des bilans de masses différentielles tout le long de l'extraction, avec les hypothèses suivants :
Cela nous permet d'obtenir les équations suivantes :
$(1-e).V.uV*\frac{\partial c_{q}}{\partial t}= -AK(q-q*)$
Les conditions initiales sont les suivantes : C = 0, q=q0 à t = 0 et c(0,t) à h=0
La phase d'équilibre est : $c = k.q*$
Sachant que le fluide et la phase sont uniformes à chaque stage, nous pouvons définir le modèle en utilisant les équations différentielles ordinaires (2n). Les équations sont les suivantes :
In [13]:
import numpy as np
from scipy import integrate
from matplotlib.pylab import *
In [14]:
import numpy as np
from scipy import integrate
import matplotlib.pyplot as plt
def vdp1(t, y):
return np.array([y[1], (1 - y[0]**2)*y[1] - y[0]])
t0, t1 = 0, 20 # start and end
t = np.linspace(t0, t1, 100) # the points of evaluation of solution
y0 = [2, 0] # initial value
y = np.zeros((len(t), len(y0))) # array for solution
y[0, :] = y0
r = integrate.ode(vdp1).set_integrator("dopri5") # choice of method
r.set_initial_value(y0, t0) # initial values
for i in range(1, t.size):
y[i, :] = r.integrate(t[i]) # get one more value, add it to the array
if not r.successful():
raise RuntimeError("Could not integrate")
plt.plot(t, y)
plt.show()
In [15]:
P = 9 #MPa
T = 323 # K
Q = 8.83 #g/min
e = 0.4
rho = 285 #kg/m3
miu = 2.31e-5 # Pa*s
dp = 0.75e-3 # m
Dl = 0.24e-5 #m2/s
De = 8.48e-12 # m2/s
Di = 6e-13
u = 0.455e-3 #m/s
kf = 1.91e-5 #m/s
de = 0.06 # m
W = 0.160 # kg
kp = 0.2
r = 0.31 #m
n = 10
V = 12
#C = kp * qE
C = 0.1
qE = C / kp
Cn = 0.05
Cm = 0.02
t = np.linspace(0,10, 1)
ti = (r ** 2) / (15 * Di)
In [16]:
def reverchon(x,t):
#Ecuaciones diferenciales del modelo Reverchon
#dCdt = - (n/(e * V)) * (W * (Cn - Cm) / rho + (1 - e) * V * dqdt)
#dqdt = - (1 / ti) * (q - qE)
q = x[0]
C = x[1]
qE = C / kp
dqdt = - (1 / ti) * (q - qE)
dCdt = - (n/(e * V)) * (W * (C - Cm) / rho + (1 - e) * V * dqdt)
return [dqdt, dCdt]
In [17]:
reverchon([1, 2], 0)
Out[17]:
In [18]:
x0 = [0, 0]
t = np.linspace(0, 3000, 500)
resultado = odeint(reverchon, x0, t)
qR = resultado[:, 0]
CR = resultado[:, 1]
plt.plot(t, CR)
plt.title("Modelo Reverchon")
plt.xlabel("t [=] min")
plt.ylabel("C [=] $kg/m^3$")
Out[18]:
In [19]:
x0 = [0, 0]
t = np.linspace(0, 3000, 500)
resultado = odeint(reverchon, x0, t)
qR = resultado[:, 0]
CR = resultado[:, 1]
plt.plot(t, qR)
plt.title("Modelo Reverchon")
plt.xlabel("t [=] min")
plt.ylabel("C solid–fluid interface [=] $kg/m^3$")
Out[19]:
In [20]:
print(CR)
In [21]:
r = 0.31 #m
x0 = [0, 0]
t = np.linspace(0, 3000, 500)
resultado = odeint(reverchon, x0, t)
qR = resultado[:, 0]
CR = resultado[:, 1]
plt.plot(t, CR)
plt.title("Modelo Reverchon")
plt.xlabel("t [=] min")
plt.ylabel("C [=] $kg/m^3$")
Out[21]:
In [22]:
r = 0.231 #m
x0 = [0, 0]
t = np.linspace(0, 3000, 500)
resultado = odeint(reverchon, x0, t)
qR = resultado[:, 0]
CR = resultado[:, 1]
plt.plot(t, CR)
plt.title("Modelo Reverchon")
plt.xlabel("t [=] min")
plt.ylabel("C [=] $kg/m^3$")
Out[22]:
In [23]:
fig,axes=plt.subplots(2,2)
axes[0,0].plot(t,CR)
axes[1,0].plot(t,qR)
Out[23]:
[1] E. Reverchon, Mathematical modelling of supercritical extraction of sage oil, AIChE J. 42 (1996) 1765–1771. https://onlinelibrary.wiley.com/doi/pdf/10.1002/aic.690420627
[2] Amit Rai, Kumargaurao D.Punase, Bikash Mohanty, Ravindra Bhargava, Evaluation of models for supercritical fluid extraction, International Journal of Heat and Mass Transfer Volume 72, May 2014, Pages 274-287. https://www.sciencedirect.com/science/article/pii/S0017931014000398
Explicaciones :
In [24]:
#Datos experimentales
x_data = np.linspace(0,9,10)
y_data = np.array([0.000,0.416,0.489,0.595,0.506,0.493,0.458,0.394,0.335,0.309])
def f(y, t, k):
""" sistema de ecuaciones diferenciales ordinarias """
return (-k[0]*y[0], k[0]*y[0]-k[1]*y[1], k[1]*y[1])
def my_ls_func(x,teta):
f2 = lambda y, t: f(y, t, teta)
# calcular el valor de la ecuación diferencial en cada punto
r = integrate.odeint(f2, y0, x)
return r[:,1]
def f_resid(p):
# definir la función de minimos cuadrados para cada valor de y"""
return y_data - my_ls_func(x_data,p)
#resolver el problema de optimización
guess = [0.2, 0.3] #valores inicales para los parámetros
y0 = [1,0,0] #valores inciales para el sistema de ODEs
(c, kvg) = optimize.leastsq(f_resid, guess) #get params
print("parameter values are ",c)
# interpolar los valores de las ODEs usando splines
xeval = np.linspace(min(x_data), max(x_data),30)
gls = interpolate.UnivariateSpline(xeval, my_ls_func(xeval,c), k=3, s=0)
xeval = np.linspace(min(x_data), max(x_data), 200)
#Gráficar los resultados
pp.plot(x_data, y_data,'.r',xeval,gls(xeval),'-b')
pp.xlabel('t [=] min',{"fontsize":16})
pp.ylabel("C",{"fontsize":16})
pp.legend(('Datos','Modelo'),loc=0)
pp.show()
In [25]:
f_resid(guess)
Out[25]:
In [26]:
#Datos experimentales
x_data = np.linspace(0,9,10)
y_data = np.array([0.000,0.416,0.489,0.595,0.506,0.493,0.458,0.394,0.335,0.309])
print(y_data)
# def f(y, t, k):
# """ sistema de ecuaciones diferenciales ordinarias """
# return (-k[0]*y[0], k[0]*y[0]-k[1]*y[1], k[1]*y[1])
def reverchon(x,t,Di):
#Ecuaciones diferenciales del modelo Reverchon
#dCdt = - (n/(e * V)) * (W * (Cn - Cm) / rho + (1 - e) * V * dqdt)
#dqdt = - (1 / ti) * (q - qE)
q = x[0]
C = x[1]
qE = C / kp
ti = (r**2) / (15 * Di)
dqdt = - (1 / ti) * (q - qE)
dCdt = - (n/(e * V)) * (W * (C - Cm) / rho + (1 - e) * V * dqdt)
return [dqdt, dCdt]
def my_ls_func(x,teta):
f2 = lambda y, t: reverchon(y, t, teta)
# calcular el valor de la ecuación diferencial en cada punto
rr = integrate.odeint(f2, y0, x)
print(f2)
return rr[:,1]
def f_resid(p):
# definir la función de minimos cuadrados para cada valor de y"""
return y_data - my_ls_func(p,x_data)
#resolver el problema de optimización
guess = np.array([0.2]) #valores inicales para los parámetros
y0 = [0,0] #valores inciales para el sistema de ODEs
(c, kvg) = optimize.leastsq(f_resid, guess) #get params
print("parameter values are ",c)
# interpolar los valores de las ODEs usando splines
xeval = np.linspace(min(x_data), max(x_data),30)
gls = interpolate.UnivariateSpline(xeval, my_ls_func(xeval,c), k=3, s=0)
xeval = np.linspace(min(x_data), max(x_data), 200)
#Gráficar los resultados
pp.plot(x_data, y_data,'.r',xeval,gls(xeval),'-b')
pp.xlabel('t [=] min',{"fontsize":16})
pp.ylabel("C",{"fontsize":16})
pp.legend(('Datos','Modelo'),loc=0)
pp.show()
In [ ]:
def my_ls_func(x,teta):
f2 = lambda y, t: reverchon(y, t, teta)
# calcular el valor de la ecuación diferencial en cada punto
r = integrate.odeint(f2, y0, x)
print(f2)
return r[:,1]
In [27]:
my_ls_func(y0,guess)
Out[27]:
In [28]:
f_resid(guess)
Out[28]:
In [ ]:
In [ ]: