#Saponification de l'éthanoate d'éthyle

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

##Etude en réacteur semi -fermé isotherme
#données
CA= 0.2     #concentration de la solution d'acétate d'éthyle  en mol/L
VA= 30      #volume de la solution d'acétate d'éthyle en L
CB= 0.3     # concentration de la solution de soude  en mol/L
Qe= 2.5     #débit volumique en entrée de la solution de soude en L/min
k=7.0       #constante de vitesse  en L/mol/min

#taux de conversion en fonction du temps t <10min

n=200       #nombre de valeurs de t
T=np.linspace ( 0,10 ,n)
def V(t):
    return VA+Qe*t
def F(X,t):
    return k*(1-X)*(CB*Qe*t-CA*VA*X)/(V(t))
Xsf = odeint ( F,[0],T)

plt.figure(0)
plt.title("taux de conversion en reacteur semi-fermé")
plt.plot(T,Xsf,'ob',markersize=1,label='Xsf(t)')
plt.legend()
plt.xlabel("t(min)")
plt.ylabel ("X")
plt.ylim(0,1)
plt.show()

Xsf10=Xsf[-1]
print ("la valeur du taux de conversion au bout de 10 min est  ",Xsf10)


# taux de conversion pour t > 10minutes
Vtot=V(10)
nA=CA*VA                #quantité totale de A introduite
nB=CB*Qe*10             #quantité totale de B introduite


def F2(x,t):
    return k*(1-x)*(nB-nA*x)/Vtot
Xf=odeint(F2,Xsf10,T)


plt.figure(1)
plt.title("taux de conversion ")
plt.plot(T,Xsf,'ob',markersize=1,label='Xsf(t)')
plt.plot(T+10,Xf,'xg',markersize=1,label='Xf(t)')
plt.plot([0,20],[0.95,0.95],'r',label='XA=95%')
plt.legend()
plt.xlabel("t(min)")
plt.ylabel ("XA")
plt.grid ()
plt.ylim(0,1)
plt.show()

#détermination de la durée pour laquelle XA=0,95
i=0
while Xf[i] <= 0.95:
    i=i+1
print ("pour avoir un taux de conversion de ", Xf[i][0]*100 ," %,la durée est de ",10+T[i],"minutes" )



##Etude en réacteur fermé adiabatique
DrH= -48657    #enthalpie standard de réaction en J/mol
Cp = 4180       #capacité thermique du mélange  en J/kg
mu= 1           #masse volumique de l'eau en kg/L
Ea = 103818
A = 1.75e19
R=8.314

#conditions  initiales
CA0=nA/Vtot
CB0=nB/Vtot
Te = 21.6+273

#Coefficient adiabatique de température
J = -DrH*CA0/Cp/mu

tad= np.linspace (0, 20 ,2*n)

# Expression de la constante de vitesse -Loi d'Arrhenius
def k(T):
    return A*np.exp(-Ea/R/T)


# Système d'équations différentielles
def F3(X,t):
    Xad,T =X[0],X[1]
    return (k(T)*(1-Xad)*(CB0-CA0*Xad),J*k(T)*(1-Xad)*(CB0-CA0*Xad))
sol=odeint(F3,[0,Te],tad)

Xfad=sol[:,0]
Tfad=sol[:,1]

plt.figure(2)
plt.title("réacteur fermé adiabatique ")
plt.plot(tad,Xfad,'ob',markersize=1,label='Xfad(t)')
plt.legend()
plt.xlabel("t(min)")
plt.ylabel ("X adiabatique")
plt.grid ()
plt.ylim(0,1)
plt.show()

plt.figure(3)
plt.title("réacteur fermé adiabatique-Evolution de la température  ")
plt.plot(tad,Tfad,'ob',markersize=1,label='Tfad(t)')
plt.legend()
plt.xlabel("t(min)")
plt.ylabel ("T")
plt.grid ()
plt.show()

# détermination du temps pour obtenir un taux de conversion de 95 %
i=0
while Xfad[i] <= 0.95:
    i=i+1
print ("pour avoir un taux de conversion de",Xfad[i],"% la durée est de ",tad[i],"minutes" )