# CALCUL DE L'INFLUENCE DE P ET T SUR L'ÉQUILIBRE DE LA SYNTHESE DE NH3

import numpy as np
import matplotlib.pyplot as plt
from scipy import optimize as sp # bibliotheque de la fonct. bisect (dichotomie)

# définition de la fonction K(T) constante de réaction
# deltarH° = -92,2 kJ/mol, deltarS° = -198,7 J/K/mol
# et approximation d'Ellingham
def K(T):
    return np.exp(-(-92.2e3-T*(-198.7))/(8.314*T))

# définition du tableau des pressions, températures
# du tableau à deux entrées rendement eij(T,p)
pj=np.arange(0,500,5)
pj[0]=1                  # pmin=1 bar, 0 fait diverger Q
Ti=np.arange(300,800,5)
eij=np.zeros((100,100))

# définition de la fonction f=Q-K(T) qui s'annule à l'équilibre
# x est l'avancement et aussi le rendement ici puisque n0(N2)=1 mole
def f(x):
    return (2*x)**2*(4-2*x)**2/p**2-K(T)*(1-x)*(3-3*x)**3

# définition de la fonction rendement e qui calcule par dichotomie
# la valeur de x qui annule f, donc celle de l'équilibre
def e():
    return sp.bisect(f,0,1)   # un avancement sur l'intervalle 0 à 1


# calcul du rendement-------------------------------------------------
for i in range(100):    # pour toutes les valeurs de T
    T=Ti[i]

    for j in range(100):# pour toutes les valeurs de p
        p=pj[j]
        eij[j,i]=e()    # on résout par dichotomie




# affichage du rendement en fonction de p pour différentes températures
plt.figure(figsize=(15,8))
for i in np.arange(0,100,10): # pour quelques valeurs de T
    plt.plot(pj,eij[:,i],label='T = '+str(Ti[i])+' K') # [:,i] prend toutes les
plt.legend()                                           # valeurs de eij pour
plt.grid()                                             # tous les j
plt.title("Rendement en fonction de la pression pour différentes températures en K")
plt.ylabel("Rendement")
plt.xlabel("Pression en bar")
plt.show()



# affichage du rendement en fonction de T pour différentes pressions
plt.figure(figsize=(15,8))
for j in [0,2,6,12,24,40,60,76,90]: # pour quelques valeurs de p
    plt.plot(Ti,eij[j,:],label='p = '+str(pj[j])+' bars')
plt.legend()
plt.title("Rendement en fonction de la température pour différentes pressions en bar")
plt.ylabel("Rendement")
plt.xlabel("Température en K")
plt.grid()
plt.show()



# rendement en fonction de T et p------------------------------------

# la fonction meshgrid construit deux tableaux T (temprérature)
# et e (rendement) de N X N valeurs permettant de faire le calcul
# de la pression de travail pour toutes les valeurs de e et de T
T,p = np.meshgrid(Ti,pj)

# affichage du rendement e en nuance de couleur : fonction contourf
# on affiche e en fonction de T et de p avec des couleurs différentes
plt.figure(figsize=(15,8))
ax=plt.axes()
ax.patch.set_facecolor('darkred')
cf=plt.contourf(T,p,eij,200,cmap='jet') # cm.jet couleurs de bleu à rouge

plt.colorbar(cf)                      # affichage de l'échelle sur le coté
plt.title("Synthèse de NH3, Valeur du rendement")
plt.ylabel("Pression en Bar")
plt.xlabel("Température en K")
plt.axis([300,790,1,450])    # pour tracer sur les bonnes gammes de T et P


c=plt.contour(T,p,eij,20,colors='black')# pour tracer les contours en noirs
plt.clabel(c)                     # avec les valeurs numériques sur les lignes
plt.show()





# étude de l'influence des proportions initiales-------------------------------
# fonction retournant nH2 initial en fonction  du ratio en gardant ntot = 4 moles
def nH2():
    return ratio*4/(1+ratio)

# idem pour nN2 initial
def nN2():
    return 4/(1+ratio)

# fonction Q-K(T) qui s'annule à l'équilibre
def f(x):
    return (2*x)**2*(4-2*x)**2/p**2-K(T)*(nN2()-x)*(nH2()-3*x)**3

# rendement trouvé par dichotomie
# attention l'intervalle de recherche est plus complexe : cela dépend
# du réactif limitant... nous ne sommes plus en proportion stoschiométriques !
def e():
    return sp.bisect(f,0,xmax)


p=350
T=700
# tableau des différentes valeurs de rendement et du ratio n(H2)/n(N2)
ei=np.zeros(100)
y=np.linspace(1,6,100)

plt.figure(figsize=(15,8))

for p in [100,200,300,400]: # pour quelques pressions

    for i in range(100):    # pour les différents ratios
        ratio=y[i]
        if ratio<3:         # si H2 limitant
            xmax=nH2()/3
        if ratio>=3:        # si N2 limitant
            xmax=nN2()

        ei[i]=e()

    plt.plot(y,ei,label='p ='+str(p)+' bars') # on trace
plt.legend()
plt.axvline(3,color='black',linestyle='--')
plt.title("Rendement en fonction du ratio initial n(H2)/n(N2) pour différentes pressions en bar à 700 K")
plt.ylabel("Rendement")
plt.xlabel("Ratio initial n(H2)/n(N2)")
plt.show()