#Décomposition du peroxyde de di tertiobutyle#

import numpy as np
import matplotlib.pyplot as plt


# Données  cinétiques  : réaction d'ordre 1 #
#contante de vitesse : k = k0exp(-Ea/RT)#
k0=10**15             #unité : s-1#
Ea=157000             #unité : J/mol#

#Donnée thermodynamique#
DrH0=-150000         #unité : J/mol#
R = 8.314           #unité : J/K/mol#
Cp = 2100           #capacité thermique massique moyenne du système chimique  , J/K/kg #
mu = 900            # masse volumique   , kg/m3 #

#Donnée sur le réacteur : RPAC #
VR = 5.2e-4      #volume du réacteur en m3#
#Alimentation #
CRe=4e3        #concentration molaire en entrée du réactif (peroxyde) en mol/m3#
tau = 600    #temps de passage en s#
Te1=600      #-température d'entrée  en kelvin#

#1-Variation du taux de conversion X en fonction de la température  T  : aspect cinétique"
Temp= np.linspace ( 250,600,500)

#calcul des valeurs de k pour chaque température
def k(T):
    return k0*np.exp(-Ea/R/T)

#calcul du taux de conversion pour chaque température #
def X1(T):
    Y = tau*k(T)
    return Y/ ( 1+Y)

Xcin=X1(Temp)

#tracé des varaiations du taux de conversion en fonction de la température #
plt.figure(0)
plt.plot(Temp,Xcin,'+b',markersize=1,label="aspect cinétique")
plt.xlabel("T(K)")
plt.ylabel("X cinétique")
plt.grid()
plt.legend()
plt.show()

#2-Aspect thermodynamique , fonctionnement adiabatique "
J=-CRe*DrH0/Cp/mu
def X2 (T):
    return (T-Te1)/J
Xadia =X2(Temp)

plt.figure(1)
plt.plot(Temp,Xadia,'+r',markersize=1,label="fonctionnement adiabatique")
plt.xlabel("T(K)")
plt.ylabel("X adiabatique")
plt.grid()
plt.legend()
plt.show()

plt.figure(2)
plt.plot(Temp,Xadia,'+r',markersize=1,label="aspect thermodynamique")
plt.plot(Temp,Xcin,'+b',markersize=1,label="aspect cinétique")
plt.xlabel("T(K)")
plt.ylim(0,1)
plt.ylabel("X adiabatique")
plt.grid()
plt.legend()
plt.show()

#3-Aspect thermodynamique, avec refroidissement par l'eau #
S= 314e-4  #surface d'échange  en m2#
U = 80  # coefficient d'écange    , W/m2/K
Teau=293    #température de l'eau en K#

E=S*U/(CRe*VR/tau)/DrH0

def X3(T):
    return E*(Teau-T)+(T-Te1)/J
Xechange = X3(Temp)

plt.figure(3)
plt.plot(Temp,Xechange,'+r',markersize=1,label="aspect thermodynamique")
plt.plot(Temp,Xcin,'+b',markersize=1,label="aspect cinétique")
plt.xlabel("T(K)")
plt.ylim(0,1)
plt.ylabel(" taux de conversion")
plt.grid()
plt.legend()
plt.show()