#Procédé de contact - 2023/2024
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import bisect

##Variations du taux de conversion en fonction de la température
#Données thermodynamiques
R=8.314         #constante des gaz parfaits  en JK-1mol-1
H298 = -197600       #enthalpie standard de réaction en Jmol-1
S298=  - 187.8       #entropie standard de réaction  en JK-1mol-1


#Constante d'équilibre dans le cadre de l'approximation d'Elligham
#Définir la fontion K (T) qui permettra de calculer les valeurs de K° pour une température T
def K(T):
    return np.exp(-(H298-T*S298)/R/T)
print(K(718))
print(K(298))


#Courbe K°(T)
# Introduire un tableau numpy contenant  les valeurs de T à considérer : intervalle [298K -800 K], on le notera Tx
Tx=np.linspace (500,1000,200)


#Créer le tableau des valeurs de K(T) , on le notera  Ky
Ky= K(Tx)

# représentation des variations de K° en fonction de T
plt.figure(0)
plt.plot(Tx,Ky,'+b',label="constante d'équilibre")
plt.xlabel("T(K)")
plt.ylabel ("K°(T)")
plt.grid()
plt.title("variations de K° en fontion de la température ")
plt.legend()
plt.show()

# taux de conversion
# Conditions  : r = 1,33   , n(N2) = 4 n(02)  et  P fixée à P°=1 bar
#expression du quotient reactionnel  en fontion de T et r
def ntot(x,r):
    return 5+r*(1-0.5*x)
def Q(x,r):
    return ntot(x,r)*x**2/(1-x)**2/(1-0.5*x*r)

# valeurs du taux de conversion en fonction de T
#le tableau contenant les valeurs du taux de conversion sera désigné par alpha , les valeurs sont déterminées en utilisant la fonction bisect
alpha=[]
def Eq(x,T,r):
    return Q(x,r)-K(T)
for i in range(len(Tx)):
    alpha.append(bisect(Eq,1e-6,1-1e-6,args=(Tx[i],1.33)))

plt.figure(1)
plt.plot(Tx,alpha,'+b',label="P=1bar",markersize=1)
plt.grid()
plt.legend()
plt.xlabel("T(K)")
plt.ylabel("Taux de conversion")
plt.show()



##  valeurs du taux de conversion en fonction de T sans Elligham
# Définir  3  fonctions permettant de calculer les valeurs de l'enthalpie standard , l'entropie standard et la constante d'équilibre , dédignée par K2

Cp= 16.8
def H(T):
    return H298+Cp*(T-298)
def S(T):
    return S298+Cp*np.log(T/298)
def K2(T):
    return np.exp(-((H(T)-T*S(T))/R/T))


# Compléter le script afin de tracer les variations du taux de conversion en fonction de T , il sera désigné par alpha2
alpha2=[]
def Eq2(x,T,r):
    return Q(x,r)-K2(T)
for i in range(len(Tx)):
    alpha2.append(bisect(Eq2,1e-6,1-1e-6,args=(Tx[i],1.33)))


plt.figure(2)
plt.plot(Tx,alpha,'ob',label="Avec hypothèse d'Ellingham",markersize=2)
plt.plot(Tx,alpha2,'+r',label="Sans hypothèse d'Ellingham",markersize=2)
plt.grid()
plt.legend()
plt.xlabel("T(K)")
plt.ylabel("Taux de conversion")
plt.show()

# Température pour laquelle le taux de conversion est inférieur à 99%

Y99 = 0.99 + 0*Tx
plt.figure(3)
plt.plot(Tx,alpha,'ob',label="Avec hypothèse d'Ellingham",markersize=2)
plt.plot(Tx,alpha2,'+r',label="Sans hypothèse d'Ellingham",markersize=2)
plt.plot(Tx,Y99,'k',label= 'taux de conversion =99%',markersize=2)
plt.grid()
plt.legend()
plt.xlabel("T(K)")
plt.ylabel("Taux de conversion")
plt.show()


## Influence de la pression   sans approximation d 'Ellingham
def Q2(x,r,P):
    return ntot(x,r)*x**2/(1-x)**2/(1-0.5*x*r)/P
alphaP1=[]
alphaP2=[]
alphaP3=[]
def Eq3(x,T,r,P):
    return Q2(x,r,P)-K2(T)
for i in range(len(Tx)):
    alphaP1.append(bisect(Eq3,1e-6,1-1e-6,args=(Tx[i],1.33,1)))
    alphaP2.append(bisect(Eq3,1e-6,1-1e-6,args=(Tx[i],1.33,5)))
    alphaP3.append(bisect(Eq3,1e-6,1-1e-6,args=(Tx[i],1.33,10)))

plt.figure(4)
plt.plot(Tx,alphaP1,'ob',label="P=1bar",markersize=1)
plt.plot(Tx,alphaP2,'+r',label="P=5bar",markersize=1)
plt.plot(Tx,alphaP3,'+g',label="P=10bar",markersize=1)
plt.grid()
plt.legend()
plt.xlabel("T(K)")
plt.ylabel("Taux de conversion")
plt.title("Influence de la pression, r = 1,33")
plt.show()


## Influence de la  composition initiale  , pression fixée à P= 1 bar
T2=np.linspace (600,1000,500)
alphar1=[]
alphar2=[]
alphar3=[]

for i in range(len(T2)):
    alphar1.append(bisect(Eq2,1e-6,1-1e-6,args=(T2[i],1.33)))
    alphar2.append(bisect(Eq2,1e-6,1-1e-6,args=(T2[i],2)))
    alphar3.append(bisect(Eq2,1e-6,1-1e-6,args=(T2[i],1)))

plt.figure(5)
plt.plot(T2,alphar1,'ob',label="r= 1.33",markersize=1)
plt.plot(T2,alphar2,'+r',label="r=2",markersize=1)
plt.plot(T2,alphar3,'xg',label="r= 1",markersize=1)
plt.grid()
plt.legend()
plt.xlabel("T(K)")
plt.ylabel("Taux de conversion")
plt.title("Influence de la composition initiale, P=1bar")
plt.show()