Accueil

 Comment faire...
avec Sa

 Tous les exercices
Sa

 Haut de page

Version PDF (69 ko)

Retour au cours III.2.2

 Retour

Sa - exercice 7

Réaction bimoléculaire réversible

A + B   ⇄  C + D      k₀, k₁     (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, r_norm

- 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 x_e. 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 A_e, 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.

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

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 A₀, les autres concentrations initiales étant fixées, ainsi que les constantes de vitesse.

Nous allons donc devoir affecter A₀ à 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 B₀, C₀, D₀ et la concentration initiale courante A₀ 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 k₀  et k₁, et nous utiliserons p[2], p[3] et p[4] pour les concentrations initiales B₀, C₀ et D₀, 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.

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 A₀

ca1[1à3][0] reçoivent p[2à4], c'est-à-dire B₀, C₀ et D₀.

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.

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 :

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] à B₀, 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 C₀ et D₀ identiques et serez ainsi ramené au cas décrit dans le cours. Vous obtiendrez deux droites en simulation.