mise à jour :31-03-2010 cinet.chim

III.2.2  Réaction bimoléculaire réversible

logo Sa

Exercice 7


Sa - exercice 7

Réaction bimoléculaire réversible

A + B   ⇄  C + D      k0, k1     (r 10)


1. Etude cinétique
bimol_rev1.cpp    bimol_rev2.cpp
 
2. Etude à l'équilibre
bimol_rev_eq.cpp

1. Etude cinétique


Nous nous proposons de calculer au cours du temps :

- les concentrations de A, B, C et D, les concentrations initiales étant quelconques

- la vitesse normalisée par la vitesse initiale, rnorm

- l'écart à l'équilibre, E

- l'avancement normalisé, χ

Jusqu'ici, nous avions toujours le choix entre la solution algébrique, en utilisant la forme intégrée des équations cinétiques, ou l'intégration numérique, avec la ou les équations différentielles, sans pouvoir préférer vraiment l'une ou l'autre.  Ici, l'utilisation de la forme intégrée (52) pose un problème : il faudrait préalablement calculer les concentrations à l'équilibre afin de déterminer p et xe. Cela serait tout à fait possible puisqu'il ne s'agit que de la solution d'une équation du second degré. Toutefois, les expressions à manipuler sont assez lourdes et les risques d'erreur d'écriture deviennent dès lors importants. De plus, les expressions compliquées sont peu lisibles et rendent la détection d'éventuelles erreurs difficile. Devant la simplicité d'écriture des équations différentielles, il est clair que la méthode à adopter dans ce cas est l'intégration numérique, comme cela le sera d'ailleurs presqu'exclusivement par la suite.

Reste à décider si on utilise une seule équation différentielle ou les quatre. Du point de vue du calcul, il serait plus logique d'en utiliser une seule, puisque le problème est monovariable. Mais, toujours pour des raisons de simplicité d'écriture, on peut préférer utiliser les quatre, sachant que quatre variables c'est encore peu pour l'intégrateur numérique. Nous proposons cette solution (télécharger bimol_rev1.cpp), que vous pouvez comparer à la solution à une équation différentielle (télécharger bimol_rev2.cpp). Le fichier de commande est identique pour les deux versions (télécharger bimol_rev.sac).

 La programmation ne fait appel qu'à des choses déjà vues. Nous attirons cependant l'attention sur le fait que le calcul de l'écart à l'équilibre, E, et de l'avancement normalisé, χ, n'est correct que si ca[0][npt-1] est suffisamment proche de Ae, c'est-à-dire si l'on peut considérer que la réaction est terminée :
for (int i = 0; i < npt; ++i)
{    
     ...
    ca[5][i] = ca[0][i] - ca[0][npt-1];  //  E
    ca[6][i] = (ca[0][0] - ca[0][i]) / (ca[0][0] - ca[0][npt-1]);  // chi

Il faudra donc veiller à ce que ce soit le cas, en attribuant à Fin une valeur suffisante,  et en vérifiant dans le fichier .sad que les concentrations ne varient pratiquement plus, pour chaque jeu de paramètres et concentrations initiales. En fait, c'est le prix à payer pour la simplicité d'écriture. Cela n'affecte en rien la précision du calcul des concentrations, qui est le principal.

Cette inexactitude du calcul de E peut fausser sa valeur relative au voisinage de zéro et donner dans ce cas des tracés en  échelle Ylog n'ayant rien à voir avec la réalité.

Inspirez-vous si nécessaire des paramètres utilisés sur la Fig. III.16 et vérifiez par simulation les points suivants :  

A- Invariance de type 2 : en divisant par deux, par exemple, toutes les concentrations initiales

     - l'échelle de temps doit être multipliée par le même facteur

     - la position de l'équilibre doit être divisée par le même facteur

     - moyennant un changement d'échelle adéquat, les courbes sont superposables.


B- Vérifier que la relaxation est une exponentielle décroissante lorsque k1 = k0. Utilisez pour cela la courbe de l'écart à l'équilibre, E, en fonction du temps. Mettez-vous dans une situation où E est positif et tracez-le en échelle Y logarithmique. Si E est négatif, vous pouvez toujours utiliser la méthode des t1/2.


C- Avec k1≠ k0, faites varier les concentrations initiales de façon à changer le signe de E. Observez que la cinétique tend à "traîner" davantage quand E est négatif. Comparez les courbes vitesse normalisée / avancement normalisé.


2. Etude à l'équilibre

télécharger bimol_rev_eq.cpp     télécharger bimol_rev_eq.sac

Nous allons maintenant voir une nouvelle possibilité très utile du logiciel Sa : l'intégration numérique automatique jusqu'à l'équilibre, ou plus exactement l'état stationnaire, qui se confond avec l'équilibre dans le cas présent. La programmation en est un peu plus difficile. Si celle-ci vous déroute un peu, vous pourrez toujours y revenir plus tard. Attachez-vous surtout à en comprendre le principe.

L'idée est d'obtenir directement la position de l'équilibre en fonction d'un paramètre donné, qui sera donc la variable indépendante du problème : ce ne sera plus le temps. Nous choisirons comme variable indépendante, dans notre exemple, la concentration initiale A0, les autres concentrations initiales étant fixées, ainsi que les constantes de vitesse.

Nous allons donc devoir affecter A0 à la variable indépendante, ind[ ], et la faire varier d'une certaine valeur Début à une valeur Fin, avec un certain Incrément, et pour chaque point injecter les concentrations initiales fixes B0, C0, D0 et la concentration initiale courante A0 dans une version spéciale de l'intégrateur numérique, srkstaComplément : intégration numérique, qui intègre automatiquement jusqu'à l'état stationnaire, et ne délivre que les valeurs de cet état (la cinétique est oubliée). Nous utiliserons ca[ ][ ] pour récupérer ces valeurs à l'équilibre. Nous garderons p[0] et p[1] pour les constantes de vitesse k0  et k1, et nous utiliserons p[2], p[3] et p[4] pour les concentrations initiales B0, C0 et D0, considérées comme paramètres du problème.

La fonction eqdiff est rigoureusement identique à celle de la version cinétique. Nous avons choisi la version à 4 équations différentielles de bimol_rev1.cpp. Les variables n_diff et first_var devront donc toujours être mises à 4 et 0, respectivement.

Pour fappel, il y a dans Sa Editor un Template : call srksta. Il ne s'agit en fait que d'un exemple, qui doit être adapté en fonction du problème traité.

Nous donnons ci-dessous les éléments importants du  fichier bimol_rev_eq.cpp afin de le commenter de façon détaillée :


//-----------------------------------------------------
void Identification(Modele& modele)
{
   ...
   //la var. ind. n'est pas la var. d'intégration :
   integvar_is_indep = false;  
}
//----------------------------------------------------
void eqdiff (Sa_data x, Sa_data* y, Sa_data* dy)
{
   dy[0] = - p[0]*y[0]*y[1] + p[1]*y[2]*y[3]; // dA/dt
   dy[1] =   dy[0];                           // dB/dt
   dy[2] = - dy[0];                           // dC/dt
   dy[3] = - dy[0];                           // dD/dt
}
//---------------------------------------------------------
void fappel()
{
   // déclarations spécifiques pour srksta  :
   Sa_data ca1[NVAR][3];
   Sa_data ind1[3];


   // boucle sur l'ensemble des points à calculer
   for (int i = 0; i < npt; i++)
   {
      // re-initialisation pour chaque point :
      for (int j = 0; j < n_diff; ++j)
         for (int m = 0; m < 3; ++m)
            ca1[j][m] = 0;

      ca1[0][0] = ind[i];  // affectation de ind[i] à A0
      ca1[1][0] = p[2];    // B0
      ca1[2][0] = p[3];    // C0
      ca1[3][0] = p[4];    // D0

      // appel de srksta. L'état stationnaire est retourné dans ca1[i][2] :
      srksta(n_diff, ca1, ind1, 0, 0, 0, h0, tol, iset, jacob, h_compt, c_min);

      // recopie de l'état stationnaire dans la variable globale ca[][]  
      ca[0][i] = ca1[0][2];
      ca[1][i] = ca1[1][2];
      ca[2][i] = ca1[2][2];
      ca[3][i] = ca1[3][2];
   }
}

Dans la fonction Identification, nous avons affecté false à integvar_is_indep. Cela indique à Sa que la variable indépendante n'est pas la variable d'intégration numérique, et que, par conséquent, il n'y a aucun lien entre l'Incrément de la variable indépendante et le pas max d'intégration. Cela rend impossible, notamment, le mode pas max Auto. Nous pourrions faire la même chose de manière interactive à l'exécution, en désactivant la case correspondante (onglet Général, volet Var. indépendante, juste sous le nom de celle-ci) , mais il est bien préférable de le prévoir dans le programme lui-même.

Toutes les variables globales de Sa, comme n_diff, firs_var, ca, ind, p, integvar_is_indep, etc. sont par définition accessibles dans votre fichier MOD.cpp. Leur déclaration, et par conséquent leur liste complète est dans le fichier ...\Sa3\Sa302_Fr\Sources\global.h (qu'il ne faut pas modifier sans savoir précisément ce qu'on fait). Toutes ont des valeurs par défaut, correspondant aux applications les plus courantes.

Le choix de l'endroit où on les initialise dépend de ce que l'on veut faire, la seule règle absolue étant de les initialiser avant de s'en servir. Ainsi, n_diff ou integvar_is_indep, par exemple, pourraient aussi bien être initialisées dans fappel : elles seraient alors ré-initialisées à chaque simulation. En les initialisant dans Identification, elles sont initialisées uniquement au lancement de l'exécutable Sa, ce qui est suffisant dans le cas présent.

fappel commence par la déclaration de deux variables locales, de type Sa_data., c'est-à-dire le type numérique par défaut de Sa. Nous les avons nommées ca1 et ind1, parce qu'elles vont jouer vis à vis de l'intégrateur srksta un rôle analogue à celui de ca et ind dans l'intégration "normale" avec srkvi. Comme ca, ca1 est un tableau à deux dimensions de NVAR lignes et 3 colonnes. NVAR (en majuscules) est le nombre de lignes effectif de ca (205 par défaut). Vous n'aviez pas à vous en servir jusqu'ici, ca étant déclaré globalement pour tout Sa. Mais ici, il est impératif de déclarer ca1 avec le même nombre de lignes. Quant au nombre de colonnes, qui correspondent à la variable d'intégration, 3 suffisent, comme pour la variable ind1. En effet, srksta va utiliser ces variables comme des variables "glissantes", ne retenant que les 3 dernières valeurs calculées. C'est d'ailleurs lorsqu'il trouvera que ces 3 valeurs sont égales, à la tolérance près, qu'il décidera d'arrêter l'intégration et de délivrer l'état stationnaire. En fait, ind1 ne sera pas réellement utilisée ici, mais elle doit être déclarée et transmise comme argument à srksta.

Il faut ensuite faire une boucle générale sur le nombre de points à calculer. L'intégrateur srksta sera appelé pour chacun des points, au cœur de cette boucle. Il faut donc tout d'abord initialiser les variables qu'il utilise.

La double boucle

   for (int j = 0; j < n_diff; ++j)
      for (int m = 0; m < 3; ++m)
         ca1[j][m] = 0;

est une pré-initialisation générale qui met à zéro tous les éléments utilisés de ca1. C'est très important car, sinon, on enverrait au départ à srksta des valeurs venant d'un calcul précédent, qui perturberaient sa détection de l'état stationnaire. C'est d'ailleurs une règle d'or en programmation que de toujours initialiser les variables à une valeur connue, et donc contrôlable.

ca1 étant ainsi "nettoyée", il faut maintenant affecter à sa 1ère colonne les vraies valeurs initiales de l'intégration :

ca1[0][0] reçoit ind[i], la valeur courante de la variable indépendante A0

ca1[1à3][0] reçoivent p[2à4], c'est-à-dire B0, C0 et D0.

Et srksta peut donc être appelé. Sa liste d'arguments ressemble à celle de srkvi, avec quelques différences à respecter à la lettre.

Il ne reste plus qu'à récupérer les valeurs des états stationnaires dans la variable globale ca. Noter que l'on pourrait recopier n'importe quelle colonne de ca1, puisqu'elles doivent être égales, mais il est tout de même plus logique de prendre la dernière.

Simulez ce système avec tous les paramètres du fichier bimol_rev_eq.sac, en mettant l'Incrément en mode log, pour être en accord avec l'Incrément, qui est de 1.1 dans le fichier, ou en changeant de façon adéquate ce dernier.

Vous aurez un message d'erreur si vous oubliez de faire comme indiqué ci-dessus. Corrigez et relancez la simulation.

Le mode log a été choisi ici parce qu'on fait varier A0 de plusieurs décades.

La première chose à faire est de vérifier que la tolérance d'intégration, qui sert aussi à la détection de l'état stationnaire, et le pas max sont corrects. Pour cela, divisez par 10, par exemple, la tolérance : vous constatez que le résultat est le même, et pouvez donc remettre la valeur initiale. Faites de même avec pas max : vous constatez que la valeur proposée était correcte et pouvez y revenir également. Nous insistons sur la nécessité de tels tests : il serait absurde de discuter sur des courbes numériquement fausses.

Ceci étant fait, vous devez observer des courbes semblables à la figure 1 :

simulation bimol_rev_eq

Fig. 1

Déplacement de l'équilibre en fonction de A0.

Résultat de simulation avec bimol_rev_eq.cpp et les paramètres de bimol_rev_eq.sac :

k0 = 1 mol−1.L.s−1 ;  k1 = 0.2 mol−1.L.s−1 

B0 = 3×10−3 mol.L−1 ;  C0 = 4×10−3 mol.L−1 ;  D0 = 3.5×10−3 mol.L−1 

Nous avons ajouté sur cette figure les horizontales représentant les concentrations initiales B0, C0 et D0 constantes, la bissectrice représentant la variation de A0, et la verticale au point d'abscisse E0 correspondant à la situation où l'écart initial à l'équilibre est nul.

- La première remarque est que la position de l'équilibre n'est pas proportionnelle à la quantité de matière mise en jeu, dans le cas général. Si c'était le cas, les courbes se réduiraient à un faisceau de droites.

- Les courbes Ce et De sont parallèles, tandis que la courbe Be est leur image retournée. C'est la traduction du bilan matière du système.

- Lorsque A0 devient grand, la courbe Ae tend vers une droite parallèle à la bissectrice : A n'est pratiquement pas consommé par la réaction parce que en grand excès.

- La courbe Be décroît quand A0 (et Ae) croît : en augmentant la vitesse directe, r0 = k0AB, la présence de A tend à "chasser" son co-réactif B.

- La verticale en E0 partage le plan en deux régions :

     (a) A0 < E0 : pour atteindre l'équilibre, A et B doivent augmenter tandis que C et D doivent diminuer. Cela correspond à un écart à l'équilibre, E, négatif (avec les conventions du cours), c'est-à-dire une évolution dans le sens inverse de la plus grande flèche (k0 > k1).

     (b) A0 > E0 : pour atteindre l'équilibre, A et B doivent diminuer tandis que C et D doivent augmenter. Cela correspond à un écart à l'équilibre, E, positif, une évolution dans le sens de la plus grande flèche.

Essayez d'autres combinaisons de paramètres. Vous ne pouvez pas simuler le cas où tous les réactifs sont en proportions stœchiométriques (paragraphe du cours) avec le programme tel que présenté ci-dessus. Il suffit cependant pour cela, par exemple, d'affecter également ind[i] à B0, l'initialisation de ca1 devient alors :
   ca1[0][0] = ind[i];  // affectation de ind[i] à A0
   ca1[1][0] = ind[i];  //  et à B0
   ca1[2][0] = p[3];    // C0
   ca1[3][0] = p[4];    // D0

(vous pouvez utiliser le même fichier .sac, le paramètre p[2] ne sera plus utilisé mais cela n'a pas d'importance)

Vous mettrez ensuite des concentrations initiales C0 et D0 identiques et serez ainsi ramené au cas décrit dans le cours. Vous obtiendrez deux droites en simulation.