Calcul des constantes d'équilibre

2.a. Code du module

Module scilab : télécharger.

Fichier de données : télécharger. La syntaxe des indices et exposants est celle de TeX/LaTeX.

La fonction readData lit le fichier de données et renvoit les données sous forme d'un objet struct.

function data=readData(filename)
    f=mopen(filename)
    n=0
    data=struct()
    while (n<>-1),
        [n,nom,h,s]=mfscanf(f,"%s%f%f\n");
        if (n<>-1) then
            data(nom)=struct('H',h,'S',s); 
        end;
    end
    mclose(f)
endfunction    
            

La fonction somme calcule l'enthalpie standard de formation et l'entropie standard pour une somme de composés affectés de coefficients. Il faut placer un * entre le coefficient et le nom. Voici par exemple une somme : C(s)+0.5*O_2(g).

function [H,S]=somme(data,som)
  sum = strsubst(som,"+}","@");
  sum = strsubst(sum,"^+","^@");
  c=tokens(som,'+');
  n=size(c);
  n=n(1);
  H=0;
  S=0;
  for i=1:n,
      x=tokens(c(i),"*");
      nx=size(x);
      nx=nx(1);
      if nx==1 then
        coef=1;
        nom=x(1);
      else
        coef=strtod(x(1));
        nom=x(2);
      end
      nom = strsubst(nom,"@","+}");
      nom = strsubst(nom,"^@","^+");
      H=H+coef*data(nom).H;
      S=S+coef*data(nom).S;
  end
endfunction
            

La fonction reaction calcule l'enthalpie standard et l'entropie standard de réaction. L'équation de la réaction est de la forme C(s)+0.5*O_2(g)=CO(g). L'enthalpie est donnée en kJ/mol et l'entropie en J/K/mol.

function [DH,DS]=reaction(data,equation)
  n=strindex(equation,"=");
  som=strsplit(equation,[n-1 n]);
  reactifs=som(1);
  produits=som(3);
  [DH,DS]=somme(data,produits);
  [h,s]=somme(data,reactifs);
  DH=DH-h;
  DS=DS-s;
endfunction
            

La fonction constanteK calcule la constante d'équilibre pour une enthalpie standard de réaction, une entropie standard de réaction et une température fournies en argument. Le calcul effectué suppose que les grandeurs standard (enthalpie et entropie) ne dépendent pas de la température.

function K=constanteK(DH,DS,T)
  R=8.31;
  K=exp(-(DH*1e3-T*DS)/(R*T)) 
endfunction
            

La fonction plotK effectue le tracé du logarithme de la constante d'équilibre en fonction de la température :

function plotK(DH,DS,Tmin,Tmax,N)
    T=[Tmin:(Tmax-Tmin)/(N-1):Tmax];
    K=T;
    R=8.31;
    for i=1:N,
        K(i)=-(DH*1e3-T(i)*DS)/(R*T(i));
    end;
    plot2d(T,K)
    xtitle("Constante d''equilibre",'T (K)','ln(K)')
endfunction
            

2.b. Exemple

On considère comme exemple la réaction suivante :

$$2S{O}_2\left(g\right)+O_2\left(g\right)=2S{O}_3\left(g\right)$$

On commence par charger le module puis le fichier de données, en indiquant son chemin par rapport au répertoire courant de scilab (le répertoire d'où le programme est lancé).

exec('thermochim.sci')
data=readData('thermodata.txt');
            

Pour accéder par exemple à l'enthalpie standard de formation du dioxyde de soufre (à 298 K), on écrira :

data('SO_2(g)').H

-296.82999 

Les grandeurs standard de la réaction sont obtenues par :

[DH,DS]=reaction(data,'2*SO_2(g)+O_2(g)=2*SO_3(g)')

[DH,DS]

-197.78003 -188.05798 

La constante d'équilibre à 700 K :

constanteK(DH,DS,700)

86684.116 

La courbe donnant la constante d'équilibre en fonction de la température :

plotA=scf();
plotK(DH,DS,300,1500,100)

Figure pleine page

3.a. Classe Thermochimie

import re
from pylab import *
class Thermochimie:

Le constructeur lit le fichier de données et construit un dictionnaire :

    def __init__(self,dataFile):
        f = open(dataFile,"r") 
        self.data = {}
        for line in f.readlines():
            li = line.strip()
            t = re.split('[\s]+',li)
            if len(t)==3:
                self.data[t[0]]={'H':float(t[1]),'S':float(t[2])}
        f.close()

La méthode somme calcule l'enthalpie standard de formation et l'entropie standard pour une somme de composés affectés de coefficients. Il faut placer un * entre le coefficient et le nom. Voici par exemple une somme : C(s)+0.5*O_2(g).

    def somme(self,som):
        sum = som.replace("+}","@")
        sum = sum.replace("^+","^@")
        cliste = re.split("\+",sum)
        H = 0
        S = 0
        for c in cliste:
            a = re.split("\*",c)
            if len(a)==2:
                coef = float(a[0])
                nom = a[1]
            else:
                coef = 1
                nom = a[0]
            nom = nom.replace("@","+}")
            nom = nom.replace("^@","^+")
            if self.data.has_key(nom):
                H = H + coef*self.data[nom]['H']
                S = S + coef*self.data[nom]['S']
            else:
                print "Erreur : le compose %s n'est pas dans la table de donnees"%nom
                raise SystemExit
        return [H,S]

la méthode reaction calcule l'enthalpie standard et l'entropie standard de réaction. L'équation de la réaction est de la forme C(s)+0.5*O_2(g)=CO(g). L'enthalpie est donnée en kJ/mol et l'entropie en J/K/mol.

    def reaction(self,equation):
        [reactifs,produits]=re.split("=",equation)
        [Hp,Sp]=self.somme(produits)
        [Hr,Sr]=self.somme(reactifs)
        return [Hp-Hr,Sp-Sr]

Exemple d'utilisation de la classe :

from thermochimie import *
th = Thermochimie("thermodata.txt")
[DH,DS]=th.reaction("2*SO_2(g)+O_2(g)=2*SO_3(g)")
                

print([DH,DS])
--> [-197.7800000000001, -188.058]

3.b. Classe EnthalpieLibre

La classe EnthalpieLibre permet de calculer et tracer l'enthalpie libre standard de réaction en fonction de la température pour différentes réactions d'oxydation.

Les équations des réactions d'oxydation sont données par le fichier texte liste_reations.txt avec le domaine de température correspondant, un caractère (O ou N) indquant si l'équation doit être écrite et un numéro précisant la réaction d'oxydation.

class EnthalpieLibre:

Le constructeur prend en argument une instance de la classe Thermochimie et le nom du fichier comportant la liste des réactions. Un dictionnaire est constuit, dont chaque élément (clé entière) est une liste d'équations correspondant à une réaction d'oxydation.

    def __init__(self,thermo,reactionsFile):
        self.thermo = thermo
        self.reaction = {}
        f = open(reactionsFile,"r")
        for line in f.readlines():
            li = line.strip()
            t = re.split('[\s]+',li)
            if len(t)==5:
                n = int(t[3])
                equation = t[0]
                texte = t[2]
                marker = t[4]
                a = re.split('=|,',t[1]) 
                Tmin = float(a[1]) 
                Tmax = float(a[2])
                if not self.reaction.has_key(n):
                    self.reaction[n] = []
                self.reaction[n].append([equation,Tmin,Tmax,texte,marker])
        f.close()

La méthode plot trace l'enthalpie libre standard en fonction de la température pour la réaction d'oxydation dont le numéro est donné en argument :

    def plot(self,n,color):
        for e in self.reaction[n]:
            equation = e[0]
            Tmin = e[1]
            Tmax = e[2]
            texte = e[3]
            mark = e[4]
            [DH,DS] = self.thermo.reaction(equation)
            G1 = DH-DS*1e-3*Tmin
            G2 = DH-DS*1e-3*Tmax
            plot([Tmin,Tmax],[G1,G2],color=color)
            if mark=="T":
                scatter([Tmax],[G2],color=color,marker='o')
            if texte=='O' or texte=='o':
                eq = equation.replace("*","")
                text(Tmax,G2,r'$\rm %s$'%eq,color=color,size=12)

La méthode gstandard calcule l'enthalpie libre standard d'une réaction pour la température demandée :

    def gstandard(self,n,T):
        G = 0
        for e in self.reaction[n]:
            equation = e[0]
            Tmin = e[1]
            Tmax = e[2]
            if (T>=Tmin) and (T<=Tmax):
                [DH,DS] = self.thermo.reaction(equation)
                G = DH-DS*1e-3*T
        return G

Exemple :

G = EnthalpieLibre(th,"liste_reactions.txt")
            

print(G.reaction[6])
--> [['2*H_2(g)+O_2(g)=2*H_2O(l)', 300.0, 373.0, 'N', 'T'], ['2*H_2(g)+O_2(g)=2*H_2O(g)', 373.0, 2000.0, 'O', 'F']]

figure(figsize=(10,6))
xlabel('T (K)') 
ylabel(r'$\Delta_rG^0\ (kJ/mol)$')
G.plot(6,'r')
G.plot(5,'b')
G.plot(8,'g')
axis([300,3000,-1000,0])
grid(True)
            

plotB.pdf

print(G.gstandard(6,1000))
--> -394.162