#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
def K(T):
    return np.exp(-(H298-T*S298)/R/T)
print(K(718))
print(K(298))

#Courbe K°(T)
T=np.linspace (500,1000,100)
Y= K(T)
plt.figure(0)
plt.plot(T,Y,'+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     r = 1,33   , n(N2) = 4 n(02)  et  P fixée à P°=1 bar
#expression du quotient reactionnel  en fontion de Tet 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 pour  r=1.33 et P fixée à P°=1bar
alpha=[]
def Eq(x,T,r):
    return Q(x,r)-K(T)
for i in range(len(T)):
    alpha.append(bisect(Eq,1e-6,1-1e-6,args=(T[i],1.33)))

plt.figure(1)
plt.plot(T,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
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))

alpha2=[]
def Eq2(x,T,r):
    return Q(x,r)-K2(T)
for i in range(len(T)):
    alpha2.append(bisect(Eq2,1e-6,1-1e-6,args=(T[i],1.33)))

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

print ()

## 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(T)):
    alphaP1.append(bisect(Eq3,1e-6,1-1e-6,args=(T[i],1.33,1)))
    alphaP2.append(bisect(Eq3,1e-6,1-1e-6,args=(T[i],1.33,5)))
    alphaP3.append(bisect(Eq3,1e-6,1-1e-6,args=(T[i],1.33,10)))

plt.figure(3)
plt.plot(T,alphaP1,'ob',label="P=1bar",markersize=1)
plt.plot(T,alphaP2,'+r',label="P=5bar",markersize=1)
plt.plot(T,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(4)
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()