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

V.4.2  Approximation de l'équilibre rapide (AER)

logo Sa

Exercice 14


Sa - exercice 14
Approximation de l'équilibre rapide (AER)

AER exemple 1
AER_bimol_1.cpp
AER exemple 2
AER_bimol_2.cpp

Cet exercice concerne l'application de l'AER lorsque la réaction réversible est bimoléculaire (au moins dans un sens). Dans les deux exemples qui suivent, bien qu'ils soient très voisins, l'AER ne s'applique pas tout à fait de la même manière.

Il permet également d'observer la fonction de stockage que peut jouer une réaction réversible, mono ou bimoléculaire, au sein d'un mécanisme.


AER exemple 1


Soit le mécanisme, analogue à (r 7abc), mais dont la réaction réversible est bimoléculaire (contre monomoléculaire) :

A  →  B     k1     (r 1a)
B + C  ⇄  D     k2, k3     (r 1b)
B  →  E     k4     (r 1c)

on suppose que seules les concentrations initiales A0 et C0 sont différentes de zéro.

Equations

Les équations exactes sont

dA/dt  =  − k1A     (1)
   dB/dt  =  k1A − k2B.C + k3D − k4B     (2)
   dC/dt  =  − k2B.C + k3D     (3)
   dD/dt  =  k2B.C − k3D     (4)
   dE/dt  =  k4B     (5)

Si les vitesses de la réaction (r 1b) sont grandes devant les autres, on peut appliquer l'AER :

D / (B.C)  ≈  k2/k3  =  K    ∀ t     (6)

L'intermédiaire total impliqué dans la réaction réversible est
   Btotal  =  B + D     (7)
on a donc, en combinant (6) et (7)
   B  =  Btotal / (1 + K.C)     (8)
   D  =  K.C.Btotal / (1 + K.C)     (9)

On a d'autre part, C ne pouvant se trouver dans d'autres espèces que D,
   C0 − C  =  D  =  K.C.Btotal / (1 + K.C)     (10)   
d'où l'équation du second degré en C :
   K.C2 + [1 + K (Btotal − C0)] C − C0  =  0     (11)
de déterminant Δ. Le produit des racines, − C0/K, étant négatif, la racine positive, physiquement acceptable est
   C  =  {− [1 + K (Btotal − C0)] + Δ1/2}/ (2K)     (12)

Le système d'équations différentielles (1) à (5) devient donc :

dA/dt  =  − k1A     (1)
   dBtotal/dt  =  k1A − k4Btotal / (1 + K.C)     (13)
   dE/dt  =  k4Btotal / (1 + K.C)     (14)

Il est intéressant d'essayer de calculer la position de l'équilibre de cette réaction composée :
les équations stationnaires, obtenues en annulant les expressions (1) à (5), donnent immédiatement
Ae  =  Be  =  0
il en résulte (éq. (4) par exemple)
De  =  0
D'autre part les équations de conservation,
A0  =  Ae + Be + De + Ee
C0  =  Ce + De
permettent de déduire :
Ee  =  A0
Ce  =  C0

Programmation

Nous proposons d'écrire un programme qui permette d'effectuer dans une seule simulation le calcul sans et avec AER, afin d'en comparer les résultats. Pour cela, nous allons approfondir quelques points de la programmation avec Sa. L'essentiel du programme est reproduit ci-dessous avec des numéros de lignes pour en faciliter les commentaires.

Télécharger  AER_bimol_1.cpp     Télécharger AER_bimol_1.sac

L'idée de base est d'utiliser les variables 0 à 5 pour le calcul exact, et les variables suivantes, de 6 à 11, pour le calcul avec AER. Toutefois, il n'est pas possible de garder exactement le même ordre des variables dans les deux cas (A, B, C, D, E, B+D pour le calcul exact ; A, B+D, E, C, B, D pour l'AER).

Tout d'abord une variable booléenne, c'est-à-dire ne pouvant prendre que la valeur true ou false, AER, est déclarée en dehors de toute fonction (ligne 1), de façon à être accessible aussi bien depuis fappel que de eqdiff.

fappel est divisé en 2 parties : l'une, initialisant AER à false (ligne 31) pour le calcul exact, l'autre, initialisant AER à true (ligne 37), pour le calcul avec AER.

La fonction eqdiff commutera (lignes 4 ou 15) en fonction de la valeur de AER vers le calcul des équations correspondantes. ca[2][0] (lignes 7 et 8) désigne la concentration initiale C0, aussi bien pour l'AER que pour le calcul exact.

Mais avant d'appeler l'intégrateur numérique, srkvi, il faut :
- affecter la bonne valeur au nombre d'équations différentielles, n_diff (lignes 32 et 38)
- permettre le décalage correct des numéros de variables pour le calcul avec AER : par défaut, la variable first_var est égale à 0, c'est pourquoi nous ne nous en sommes jamais préoccupés jusqu'à maintenant. A la ligne 39, on lui affecte la valeur 6, cela indique à l'intégrateur qu'il doit opérer sur les variables à partir du numéro 6. Noter qu'il est nécessaire de lui affecter la valeur 0, à la ligne 33, sinon, la valeur 6 resterait en vigueur lors d'un nouveau calcul. Ces affectations de n_diff et first_var dans fappel écrasent celles qui sont effectuées au lancement du programme dans la fonction Identification.
- s'assurer d'effectuer les calculs avec les mêmes concentrations initiales A0 et C0. Cela pourrait être fait à l'exécution, bien sûr, mais comme ces variables sont dupliquées, il y a un risque d'oubli de leur affecter les mêmes valeurs, et il est plus sûr de forcer les choses par programme. C'est le but des lignes 40 et 41. Noter l'utilisation de first_var dans les indices de tableaux, comme ca[first_var][0], ca[first_var+3][0] , etc..

La fonction sqrt(x) (lignes 9 et 48) retourne la racine carrée de x. La liste de toutes les fonctions mathématiques utilisables est accessible dans Sa Editor : Templates > Math functions.

Dans le calcul exact, la somme B+D est calculée après l'intégration, pour une comparaison complète avec l'AER (lignes 35-36). Noter l'absence des signes { et } dans cette boucle : ils ne sont obligatoires que lorsque la boucle comporte plusieurs instructions. 

Dans la boucle lignes 43 à 54, ce sont les concentrations "en équilibre", C, B et D, qui sont calculées (ou recalculées) après l'intégration. 

//----------------------------------------------------------------------------
 1 bool AER;
//----------------------------------------------------------------------------

 2 void eqdiff (Sa_data x, Sa_data* y, Sa_data* dy)
 3 {
 4    if (AER) // application de l'AER
 5    { 
 6       Sa_data B,b,delta,C;
 7       b = 1 + p[1]/p[2]*(y[1] - ca[2][0]);
 8       delta = b*b + 4*p[1]/p[2]*ca[2][0];
 9       C = (- b + sqrt(delta)) / (2*p[1]/p[2]);
10       B = y[1] / (1 + p[1]/p[2]*C);

11       dy[0] = - p[0]*y[0];         // dA/dt
12       dy[1] = p[0]*y[0] - p[3]*B;  // d(B+D)/dt
13       dy[2] = p[3]*B;              // dE/dt
14    }
15    else // calcul exact
16    { 
17       Sa_data r1,r2,r3,r4;

18       r1 = p[0]*y[0];        // A -> B
19       r2 = p[1]*y[1]*y[2];   // B + C -> D   
20       r3 = p[2]*y[3];        //       <-
21       r4 = p[3]*y[1];        // B -> E

22       dy[0] = - r1;                // dA/dt
23       dy[1] = r1 - r2 + r3 - r4;   // dB/dt
24       dy[2] = - r2 + r3;           // dC/dt
25       dy[3] = r2 - r3;             // dD/dt
26       dy[4] = r4;                  // dE/dt 
27    }
28 }
//----------------------------------------------------------------------------
29 void fappel()
30 {
      // 1. traitement exact (var. 0 à 5)
31    AER = false;
32    n_diff = 5;
33    first_var = 0;
34    srkvi(n_diff, &ca[first_var], ind, npt, h0, tol, iset, jacob, h_compt, c_min);
35    for (int i = 0; i < npt; i++)
36       ca[5][i] = ca[1][i] + ca[3][i];    // B+D
 
      // 2. application AER (var. 6 à 11)
37    AER = true;
38    n_diff = 3;
39    first_var = 6;
40    ca[first_var][0] = ca[0][0];   // A0
41    ca[first_var+3][0] = ca[2][0]; // C0
42    srkvi(n_diff, &ca[first_var], ind, npt, h0, tol, iset, jacob, h_compt, c_min);
43    for (int i = 0; i < npt; i++)
44    {
45       Sa_data B,D,b,delta,C;
46       b = 1 + p[1]/p[2]*(ca[first_var+1][i] - ca[2][0]);
47       delta = b*b + 4*p[1]/p[2]*ca[2][0];
48       C = (- b + sqrt(delta)) / (2*p[1]/p[2]);
49       B = ca[first_var+1][i] / (1 + p[1]/p[2]*C);
50       D = p[1]/p[2]*C*ca[first_var+1][i] / (1 + p[1]/p[2]*C);
51       ca[first_var+3][i] = C;
52       ca[first_var+4][i] = B;
53       ca[first_var+5][i] = D;
54    }
55 }
//----------------------------------------------------------------------------

Simulations


A- Faites une simulation avec les paramètres donnés dans le fichier AER_bimol_1.sac. Ils ont été choisis de sorte que l'équilibre relaxe rapidement par rapport aux autres réactions. Vous pouvez constater, par exemple en les traçant sur deux pages graphiques distinctes, que les résultats du calcul exact et de l'AER sont identiques. Ayez la curiosité d'aller voir combien de décimales sont identiques dans les données numériques (les valeurs de A ne sont pas affectées par l'AER et sont donc rigoureusement égales).


B- Observez la forme de la courbe de C : elle décroît temporairement, puis revient finalement à sa concentration initiale.

Tracez D en fonction de B. Vous constatez qu'il s'agit presque d'une droite. En réalité, on a D / (B.C) = K, mais C ne varie que faiblement, de sorte qu'on a presque D / B = cte. 


C- Diminuez la vitesse de relaxation de l'équilibre, sans changer la constante d'équilibre. Observez à partir de quand l'AER n'est plus acceptable (pour cet exemple).


D- k2 et k3 étant à la limite de l'acceptable trouvé ci-dessus, diminuez la concentration initiale de C. Bien que cela diminue encore la vitesse de relaxation de l'équilibre, vous constatez que l'AER reste acceptable. En fait, vous constatez également que la vitesse de formation de E augmente : c'est que, C diminuant, très peu de l'intermédiaire B est stocké, temporairement, sous la forme D. Le système se rapproche de  A → B → E.

A l'inverse, si vous augmentez la concentration initiale de C, la quantité de matière stockée temporairement devient importante et la formation de E est retardée.

Il s'agit là d'un mécanisme simple et, avec un peu d'expérience, on peut prévoir et comprendre en détail ce qui va se passer. La présence de réactions réversibles bimoléculaires dans un mécanisme peut cependant rendre les choses beaucoup plus difficiles. Ce serait le cas, par exemple, si C était impliqué dans d'autres réactions.


AER exemple 2


Considérons maintenant un mécanisme identique au précédent sauf par la dernière réaction : celle-ci a maintenant lieu à partir de D, et non de B.

A  →  B     k1     (r 2a)
B + C  ⇄  D     k2, k3     (r 2b)
D  →  E     k4     (r 2c)

Equations

Les équations exactes ne posent aucun problème particulier.

Pour l'application de l'AER, les équations (6) à (9) sont toujours valables. Par contre, on n'a plus cette fois  C0 − C  =  D  mais :
   C0 − C  =  D + E     (15)
en effet, dans l'exemple précédent, la matière C était simplement en équilibre avec D, sans être transformé en E à partir de D ; c'est ce qui faisait que C n'était consommé que transitoirement, et restitué en fin de réaction. Ce n'est plus le cas maintenant : C et B sont transformés en E par la voie D, et définitivement consommés.

L'équation du second degré donnant C devient donc :

  K.C2 + [1 + K (Btotal − C0 + E)] C − C0 + E  =  0     (16)
de déterminant Δ'. Le produit des racines, − C0 + E) / K, étant négatif (E ne peut être supérieur à C0), la racine positive, physiquement acceptable est
   C  =  {− [1 + K (Btotal − C0 + E)] + Δ'1/2}/ (2K)     (17)

et le système d'équations différentielles devient :

dA/dt  =  − k1A     (18)
   dBtotal/dt  =  k1A − k4 K.C.Btotal / (1 + K.C)     (19)
   dE/dt  =  k4 K.C.Btotal / (1 + K.C)     (20)

Position de l'équilibre :
Les équations stationnaires
dA/dt  =  − k1Ae  =  0     (21)
dB/dt  =  k1Ae − k2Be.Ce + k3De  =  0     (22)
dC/dt  =  − k2Be.Ce + k3De  =  0     (23)
dD/dt  =  k2Be.Ce − k3De − k4De  =  0     (24)
dE/dt  =  k4De  =  0     (25)
donnent immédiatement
Ae  =  De  =  0
il en résulte  BeCe  =  0     (26)
Les équations de conservation
A0  =  Ae + Be + De + Ee
et (15) donnent les relations
A0  =  Be + Ee
C0  =  Ce + Ee
qui, reportées dans (26) donnent l'équation
(A0 − Ee)(C0 − Ee)  =  0     (27)
Il y a donc deux solutions :
1) Ee  =  A0  ;  Be  =  0  ;  Ce  =  C0 − A0
possible si  C0 ≥ A0
2) Ee  =  C0  ;  Ce  =  0  ;  Be  =  A0 − C0
possible si  A0 > C0

La programmation est analogue à celle de l'exemple 1.

Télécharger AER_bimol_2.cpp     Télécharger AER_bimol_2.sac


Simulations


A- Simulez avec les paramètres du fichier AER_bimol_2.sac et comparez, comme dans l'exemple précédent les résultats obtenus par le calcul exact et par l'AER.


B- Observez la forme de la courbe de C, très différente de l'exemple 1.


C- Comme pour l'exemple 1, trouvez les limites d'acceptabilité de l'AER.


D- Mettez k2 = 102 mol−1.L.s−1 et k3 = 0.1 s−1. L'AER doit être acceptable. Diminuez la concentration initiale de C, par exemple C0 = 5×10−4 mol.L−1.

L'AER est-elle toujours aussi acceptable ?

Sous quelle forme se trouve la matière en fin de réaction ?


E- Gardez les valeurs ci-dessus, sauf k3 = 2 s−1, de façon à modifier la constante d'équilibre en faveur de B. Observez le stockage sous B, temporaire pour une part et définitif pour une autre part, et le ralentissement global qui en résulte.