Module LFE (Lecture Format Endf )

Ce module lit les bases de données de sections efficaces, qui sont au format ENDF (Evaluated Nuclear Data File), et les interprète afin d'indiquer quelles informations se trouvent dans chacune des bases examinées. En effet, le format ENDF est assez compliqué à lire : les informations associées à un type de réaction donné sont dans un format qui diffère selon le type de réaction, ce qui permet d'avoir des bases de données très compactes, mais les rend difficiles à interpréter pour un lecteur humain. Ainsi, LFE lit les fichiers et les interprète ; il détermine quelles informations sont disponibles, et dans quel fichier elles se trouvent.

Module ACE (A Compact ENDF )

Ce module s'appuie sur les informations obtenues par LFE ; il a pour tâche de fournir les sections efficaces de chaque réaction, à toutes les énergies entre 0 et 20 MeV, pour chaque noyau particulier (Z, A, isomère) qu'on lui demande, à une température donnée. On peut lui préciser le nom de la base de donnée où doit être prise l'information, ou bien lui donner une liste de noms de bases de données où chercher, par ordre de priorité.

Avant de procéder à la construction des sections efficaces demandées, ACE vérifie, dans sa base générale, qu'il n'a pas déjà fait ce noyau à cette température, avec cette base de données. En effet, ce n'est pas la peine de refaire le calcul à chaque fois que le noyau est demandé, le résultat ne changera pas. C'est pourquoi ACE mémorise, dans une base générale, depuis le début de son utilisation dans notre équipe, les résultats pour tous les noyaux qui lui ont été demandés.

Ce module reçoit, du programme principal, une description de la géométrie du réacteur et la convertit en cellules élémentaires utilisables dans MCNP.

Pour définir les éléments de volume du réacteur, on fournit à CGM une description en « poupée russe » : on décrit l'enceinte englobante, les grands éléments contenus dans cette enceinte, les éléments à l'intérieur de chaque sous élément, etc. Cette description se fait au moyen de figures volumiques standards, sphère, cylindre, demie-sphère, cube, demi-cylindre, etc. dont on donne la position et les dimensions. C'est un format dans lequel il est assez aisé de décrire avec précision une géométrie même relativement complexe.

A partir de ces informations, CGM se charge de construire la description des volumes dans le format d'entrée de MCNP. Ce format n'est pas intuitif, il est fondé sur la définition des surfaces délimitant les volumes élémentaires, non sur la définition de volumes. CGM joue ainsi un rôle précieux de traduction d'une description de volumes à partir de formes conceptuelles vers un format qui convient à MCNP, dont le calcul s'organise autour des surfaces délimitant les éléments de volume du réacteur.

Composition des éléments de volume

Ce que l'on appelle la géométrie du réacteur, au sens MCNP, comprend d'autres informations que les seules forme et position des éléments de volume. Il faut lui fournir également la composition de chaque volume élémentaire, sa densité, sa température, les sections efficaces des différentes interactions possibles sur les noyaux dans le volume. Avant chaque calcul MCNP, une composition est fournie à CGM qui construira le fichier de commande destiné à MCNP.

Donner la composition à MCNP, c'est donner, pour chacun des noyaux présents dans un élément de volume :

  • Le numéro atomique de l'élément, Z
  • La masse atomique de l'élément, A
  • Le suffixe : un nombre (limité à deux chiffres) qui permet de distinguer plusieurs jeux de sections efficaces pour des noyaux de même Z et A, c'est-à-dire Isomère (ce sont des atomes de même Z et A qui correspondent à des états énergétiques différents). Les sections efficaces d'interaction peuvant varier d'un isomère à l'autre, il faut préciser le jeu de sections efficaces à utiliser selon l'isomère.
  • La température : les sections efficaces changent avec la température. Là aussi, il fautpouvoir spécifier le jeu de sections efficaces à utiliser.
  • Dans quelle proportion (pourcentage) cet élément (défini par Z, A, suffixe) est présent dans l'élément de volume
  • Le chemin des fichiers où se trouvent les sections efficaces d'interaction avec les neutrons (diffusion élastique, diffusion inélastique, (n,$ \gamma$), (n,2n), fission, etc) de l'élément.

Construction du fichier de commande pour MCNP

Outre les compositions, il faut préciser à MCNP quelles valeurs il doit calculer, ce que l'on appelle les tallys. La liste des tallys à demander est fournie à CGM par le module appelant (REM), comme le sont les compositions.

CGM génère le fichier de commande destiné à MCNP avant chaque calcul MCNP, à l'initiative de REM. Ce fichier comprend la description de toutes les cellules et, pour chaque cellule, sa composition et la liste des tallys à calculer.

Module LTM (Lecture de Tallys MCNP )

Les résultats des tallys demandés sont stockés par MCNP dans un fichier (son suffixe est m). LTM lit ce fichier, il met les résultats sous forme de tableau qu'il passe au module appelant (REM).

Module DPH (Données Physiques )

Certains des noyaux présents dans le volume sont radioactifs, ils ont souvent plusieurs voies de désintégration. L'équation différentielle d'évolution de la composition des éléments de volume fait intervenir les constantes de décroissances (inverse de la période) des diverses voies de désintégration. C'est le module DPH qui, à la demande du module REM, prépare ces données à partir de la base de données JEFF3T dans laquelle se trouve, pour chaque noyau, la liste de ses modes de décroissance et leur période.

Module IME (Intégration de la Matrice d'Evolution)

C'est ce module qui assure la partie déterministe du calcul pour chaque élément de volume. L'équation à intégrer est :

$ {\frac{{dN_{i}}}{{dt}}}$ = $ \Phi$$ \Sigma_{{i\neq j}}^{}$ < Yj $\scriptstyle \rightarrow$ i$ \sigma^{{cap}}_{{j}}$$ \varphi$ > Nj + $ \Sigma_{{i\neq j}}^{}$Y'j $\scriptstyle \rightarrow$ i$ \lambda_{{j}}^{}$Nj - $ \Phi$ < $ \sigma_{{i}}^{{abs}}$$ \varphi$ > Ni - $ \lambda_{{i}}^{}$Ni

avec :

$ \Phi$: le nombre de neutrons source par seconde

$ \varphi$ : le flux de neutrons par neutron source

Ni : le nombre de noyaux de l'élément i à l'instant t

Yj $\scriptstyle \rightarrow$ i : le rapport d'embranchement de la réaction menant du noyau j vers le noyau i

$ \Phi$ < Yj $\scriptstyle \rightarrow$ i$ \sigma_{{j}}^{{cap}}$$ \varphi$ > Nj : le taux moyen de réaction sur le noyau j conduisant au noyau i

Y'j $\scriptstyle \rightarrow$ i : le rapport d'embranchement de décroissance du noyau j conduisant au noyau i

$ \lambda_{{j}}^{}$ : la constante de décroissance du noyau j

$ \Phi$ < $ \sigma_{{i}}^{{abs}}$$ \varphi$ > Ni : le taux moyen de réaction d'absorption d'un neutron par le noyau i

$ \lambda_{{i}}^{}$ : la constante de décroissance totale du noyau i

On peut écrire cette équation sous forme matricielle, où le vecteur $ \overrightarrow{N} $ représente l'ensemble des noyaux d'un élément de volume et A la matrice des coefficients (sections efficaces, décroissances radioactives) :

$ {\frac{{d\overrightarrow{N}}}{{dt}}}$ = A$ \overrightarrow{N} $

Dans le cas particulier des réacteurs à sels fondus (RSF), on introduit un terme supplémentaire dans l'équation différentielle représentant l'évolution de la composition. L'équation devient :

$ {\frac{{d\overrightarrow{N}}}{{dt}}}$ = A$ \overrightarrow{N} $ + $ \overrightarrow{B} $

où le vecteur $ \overrightarrow{B} $ représente l'alimentation : il permet de spécifier le nombre de noyaux ajoutés ou retirés par seconde, dans les éléments de volume.

Les éléments des matrices d'évolution sont fournis par le module appelant (REM). Le vecteur $ \overrightarrow{N} $, la composition courante dans l'élément de volume, et le vecteur $ \overrightarrow{B} $ sont, eux aussi, fournis par REM.

IME se charge de l'intégration de l'équation différentielle. Cette intégration se fait par une méthode Runge Kutta du 4ième ordre (RK4) à pas variable, qui permet d'obtenir une intégration plus stable qu'avec les méthodes d'intégration linéaires.

Il s'agit donc de calculer un nouveau vecteur, $ \overrightarrow{N_{1}} $ qui représente la nouvelle composition au bout du temps d'intégration, disons $ \Delta$t, dont nous allons voir comment il est déterminé afin de conserver une précision de calcul suffisante.

  • IME fait une première intégration sur le temps $ \Delta$t, et obtient le vecteur $ \overline{{N_{1}}}$ qui correspond au vecteur à l'instant t+$ \Delta$t.
  • IME fait ensuite deux intégrations successives sur des temps d'intégration divisés par deux, obtenant $ \overline{{N'_{1}}}$ à t+ $ {\frac{{\Delta \textrm{t}}}{{2}}}$ et $ \overline{{N'_{2}}}$ à t+$ \Delta$t.
  • IME vérifie que l'écart quadratique entre les vecteurs $ \overline{{N_{1}}}$ et $ \overline{{N'_{2}}}$n'est pas supérieur à une valeur consigne fixée par le programme appelant.
     
    • Si la consigne est satisfaite, IME garde $ \overline{{N'_{2}}}$comme nouvelle composition, et agrandit le pas d'intégration pour l'intégration suivante.
 
    • Si la consigne n'est pas satisfaite, IME refait le calcul ci-dessus, mais en prenant comme pas d'intégration, $ {\frac{{\Delta \textrm{t}}}{{2}}}$, et comme demi temps $ {\frac{{\Delta \textrm{t}}}{{4}}}$. Il vérifie à nouveau si la consigne est satisfaite. Si oui, il garde la nouvelle valeur de $ \overline{{N'_{2}}}$ (qui correspond, dans ce cas, au temps t+ $ {\frac{{\Delta \textrm{t}}}{{2}}}$) mais n'augmente pas le pas d'intégration suivant : la composition change trop vite, il faut continuer à intégrer sur des petits pas.
 
    • Si la consigne n'est toujours pas satisfaite, IME refait le calcul en divisant de nouveau le pas d'intégration par deux, etc.

Ces vérifications sont faites cellule par cellule. Ainsi, certains éléments du volume, dans lesquels la composition change lentement, peuvent avoir des pas d'intégration assez grands et d'autres, dans lesquels la composition change rapidement, de petits pas d'intégration. Les différentes cellules sont synchronisées régulièrement, à l'occasion du calcul des paramètres globaux du réacteur (temps d'ajustement, de synchronisation).

Grâce aux vérifications d'écarts quadratiques sur les vecteurs de composition calculés, on obtient un pas d'intégration qui est ajusté à la dynamique du réacteur, cellule par cellule. Pendant les phases transitoires, le vecteur $ \overrightarrow{N} $ change rapidement, on a des temps d'intégration courts du fait que la consigne n'est pas (toujours) satisfaite. Lorsque le réacteur atteint l'équilibre, le vecteur $ \overrightarrow{N} $ change peu, les pas d'intégration peuvent être grands.

Remarque : sur des temps longs, le flux moyen calculé par MCNP ne peut plus être traité comme constant. La variation des < $ \sigma$$ \phi$> peut cependant être estimée (par REM), ce qui permet de continuer à suivre l'évolution des compositions sans repasser par un calcul MCNP, qui est de loin plus long à faire que ces intégrations RK4.

Ajustement

IME boucle sur des intégrations successives, en ajustant ses pas d'intégration (RK4) de façon indépendante dans chaque élément de volume. Mais on demande aussi à IME de vérifier que la puissance du réacteur et sa réactivité ne changent pas au cours du temps. Ceci suppose qu'IME suspende les intégrations au même instant dans toutes les cellules, qu'il calcule ces paramètres, qu'il vérifie que leur variation ne dépasse pas la consigne donnée. Ainsi, IME reçoit une consigne pour les paramètres qu'on lui demande de contrôler et aussi un intervalle de temps, appelé « temps d'ajustement ». IME modifie le « temps d'ajustement » selon la stabilité des paramètres qu'il contrôle.

Lorsque le temps d'ajustement est atteint dans toutes les cellules, IME calcule la puissance et la réactivité . Il vérifie que la variation de ces paramètres, depuis le dernier temps d'ajustement, est inférieure à la consigne qui lui a été fixée. Si l'un de ces paramètres varie trop, IME recommence ses calculs sur un intervalle de temps (temps d'ajustement) plus court.

Le prochain « temps d'ajustement » est fixé en fonction de la difficulté à respecter les consignes :

  • Si les consignes sont satisfaites d'emblée, et qu'elles l'avaient été au calcul précédent, le « temps d'ajustement » est augmenté pour le calcul suivant.
  • Si les consignes sont satisfaites d'emblée, mais qu'elles ne l'avaient pas été au calcul précédent (on avait réduit le « temps d'ajustement ») IME garde le même « temps d'ajustement » pour la suite.
  • Si les consignes n'étaient pas satisfaites d'emblée, IME prend comme nouveau « temps d'ajustement » le temps qui a permis de satisfaire les consignes.

Notons qu'on peut demander à IME de travailler à puissance constante. Dans ce cas, IME corrige directement les éventuelles variations de puissance en agissant sur le nombre de neutrons source, ce qui reviendrait, s'il s'agissait d'un réacteur hybride, à changer l'intensité du faisceau et, s'il s'agissait d'un réacteur critique, à changer le niveau de flux neutronique, par exemple en agissant sur les barres de contrôle.

  La puissance dépend du nombre de fissions qui se produisent, et de l'énergie libérée lors d'une fission. IME la calcule, sur l'ensemble des cellules, par :

P = $ \Phi$$ \Sigma_{{i}}^{}$Ni < $ \sigma_{{i}}^{{fis}}$$ \varphi$ > $ \varepsilon_{{i}}^{}$

avec :

$ \Phi$: le nombre de neutrons source. C'est sur ce terme qu'IME agit pour ajuster la puissance.

Ni : le nombre de noyaux de l'élément fissile i

< $ \sigma_{{i}}^{{fis}}$$ \varphi$ > : la section efficace moyenne de fission du noyau i

$ \varepsilon_{{i}}^{}$ ; l'énergie de fission du noyau i, c'est-à-dire l'énergie libérée lors de la fission d'un noyau de l'élément i.

  La réactivité est, par définition, le nombre de neutrons créés par les fissions divisé par le nombre total de neutrons absorbés, ce qui s'écrit :
k = $ {\frac{{n_{cr\acute{e}\acute{e}s}^{fis}}}{{n_{tot}^{abs}}}}$ = $ {\frac{{\Phi \Sigma _{i}N_{i}\upsilon _{i}<\sigma _{i}^{fis}\varphi >}}{{\Phi ... ...igma _{i}^{n,2n}\varphi >-2\Phi \Sigma _{i}N_{i}<\sigma _{i}^{n,3n}\varphi >}}}$ où $ \upsilon_{{i}}^{}$ est le nombre moyen de neutrons créés lors de la fission d'un noyau i

Quand les paramètres sont dans les limites fixées, IME appelle une procédure du module appelant (en jargon informatique, c'est un « callback ») qui lui rendra la main ultérieurement. Par l'intermédiaire de cette procédure, tout type d'ajustement peut être fait, qui simule une intervention extérieure sur le réacteur, et/ou sur l'alimentation, par exemple pour ajuster la réactivité. En effet, IME lui-même ne peut pas ajuster la réactivité, car cela suppose une modification de la composition des cellules, ou du vecteur d'alimentation. Ce genre d'action est réservé à un module de plus haut niveau, voire au programme principal.

Synchronisation

Un instant, dit « temps de synchronisation » peut être fixé à IME. Il est mesuré en temps absolu : le calcul sera suspendu à un instant ts qui est le temps écoulé depuis le démarrage du réacteur. Lorsque ce temps est atteint dans toutes les cellules, IME calcule, comme lors des ajustements, la puissance et la réactivité, puis appelle une procédure du module appelant qui lui rendra la main ultérieurement.

IME fonctionne ainsi de façon autonome, calculant les compositions des différents éléments de volume du réacteur au cours du temps, par intégrations successives de l'équation différentielle. Mais ceci ne continue pas indéfiniment. IME s'arrête au bout d'un temps fixé par le programme appelant, dit « temps total d'intégration » (également appelé « pas MCNP ») et lui rend la main. En effet, il arrive un moment où il faut refaire un calcul MCNP avec les compositions trouvées par IME, obtenir les nouvelles valeurs de < $ \varphi$ > et des < $ \sigma$$ \varphi$ >, puis reprendre les intégrations avec une matrice d'intégration mise à jour. Ceci se fait sous le contrôle du module appelant, (REM).

C'est le module qui contrôle l'ensemble du calcul d'évolution du réacteur, il dialogue avec les autres modules.


Construction de l'arbre d'évolution

REM connaît la composition de départ de chacune des cellules. A partir des données nucléaires - décroissances radioactives (obtenues de DPH), interactions possibles (obtenues de ACE) - REM peut déterminer, de proche en proche, quels noyaux vont intervenir au cours d'une simulation. Il construit ainsi un « arbre d'évolution » qui contient tous les noyaux présents initialement dans le réacteur, et ceux qui sont produits plus tard. Cet arbre peut contenir plusieurs milliers de noyaux, ce qui est beaucoup. On cherche à élaguer l'arbre afin de limiter la taille du système à intégrer (IME).

  • On peut éliminer les noyaux dont la période est inférieure à une valeur donnée. Cela revient à dire que le noyau se désintègre instantanément en son descendant, il ne figure dans l'arbre que par son fils.
  • On peut également en éliminer tous les noyaux dont la production suppose un nombre de réactions (profondeur de l'arbre) supérieur à une valeur donnée - la probabilité que de tels noyaux soient produits est considérée comme si faible que les éliminer du calcul n'a pas de conséquence significative.

Le seuil de période et la profondeur maximum de l'arbre sont des paramètres qui sont fournis à REM à son lancement. Au fur et à mesure que l'arbre se construit, ACE est invoqué, et les fichiers de section efficace d'interaction des différents noyaux sont complétés.


Création de xsdir

Une fois que l'arbre est complet, REM dit à ACE de construire le fichier xsdir qui contient la liste des noyaux et le chemin d'accès au fichier binaire des sections efficaces d'interaction du dit noyau. A ce point, tout est prêt pour lancer la boucle de calcul.


Organisation du calcul

REM fournit à CGM, pour chaque cellule du volume, la composition (le vecteur $ \overrightarrow{N} $ courant de la cellule), la température, la densité dans la cellule, les tallys, et il lui demande de construire le fichier d'appel à MCNP. Puis REM déclenche le calcul MCNP. Il invoque ensuite LTM pour récupérer les tallys que MCNP a calculés. Il dispose ainsi des < $ \sigma_{{i}}^{}$$ \varphi$ > et des < Yj $\scriptstyle \rightarrow$ i$ \sigma^{{cap}}_{{j}}$$ \varphi$ > dans toutes les cellules du volume. Il va pouvoir calculer les éléments de matrice d'évolution pour les passer à IME, en s'appuyant aussi sur les $ \lambda_{{i}}^{}$ et les Y'j $\scriptstyle \rightarrow$ i$ \lambda_{{j}}^{}$ qui viennent de DPH. REM doit passer à IME :

  • La matrice d'évolution A pour chaque cellule.
  • Le vecteur composition $ \overrightarrow{N} $ pour chaque cellule.
  • Le vecteur d'alimentation (dans le cas d'un RSF) $ \overrightarrow{B} $ pour chaque cellule.
  • Les consignes à respecter.
  • Le temps d'ajustement initial, et la procédure à invoquer chaque fois que le temps d'ajustement est atteint (le temps d'ajustement est variable, contrôlé par IME en fonction de la précision demandée).
  • L'instant de synchronisation et la procédure à invoquer quand il est atteint.
  • Le temps total d'intégration pour IME.

C'est REM qui déclenche un calcul IME.

Chaque fois qu'un temps d'ajustement est atteint, la procédure d'ajustement est invoquée. Elle mémorise les résultats (la puissance P, la réactivité k, la composition $ \overrightarrow{N} $ de chaque cellule et, dans le cas des réacteurs à sels fondus, l'alimentation $ \overrightarrow{B} $ de chaque cellule). La procédure peut invoquer une procédure du programme appelant (le programme principal) qui, lui-même, peut modifier les compositions, les alimentations. Finalement, la procédure rend la main à IME qui reprend son calcul, là où il en était, avec les nouveaux paramètres.

Lorsque la procédure de synchronisation est invoquée (le temps de synchronisation est atteint dans IME), les données sont mémorisées comme lors des ajustements. Des éléments de la matrice d'intégration, les compositions, les alimentations peuvent être ajustés par une procédure du programme appelant qui peut être invoquée comme dans le cas des ajustements. Un nouveau temps de synchronisation, et une nouvelle procédure de synchronisation sont passés à IME avant de lui redonner la main pour qu'il reprenne, là encore, le calcul là où il en était.

Il s'agit, avec ces « temps de synchronisation » de donner au programme principal la possibilité de faire des modifications des compositions, de la matrice d'intégration, du vecteur d'alimentation à des instants prédéterminés, simulant ainsi une intervention sur le réacteur et son mode de fonctionnement.

Il s'agit aussi, avec les « temps d'ajustement » et les « temps de synchronisation », de mémoriser régulièrement les paramètres caractérisant le fonctionnement du réacteur.

Quand le « temps total d'intégration » est atteint, IME rend la main à REM.


Contrôle du « pas MCNP »

Afin de limiter le nombre de calculs MCNP, (très longs par rapport aux calculs IME), il est tentant de faire une interpolation sur les < $ \sigma$$ \varphi$ > fournis à IME. Nous avons choisi d'utiliser une spline pour interpoler les < $ \sigma$$ \varphi$ > entre deux calculs MCNP, ce qui nous permet d'obtenir une dérivée continue autour d'un calcul MCNP.

sig_interpol.jpg
Exemple de spline utilisée pour l'interpolation entre deux calculs MCNP

Ainsi, on donne à IME non pas la valeur des < $ \sigma$$ \varphi$ >, mais leur évolution, c'est-à-dire les paramètres (a,b,c,d) de la fonction spline at3 + bt2 + ct + d qui approxime l'évolution des < $ \sigma$$ \varphi$ >. Encore faut-il connaître les deux valeurs de < $ \sigma$$ \varphi$ > entre lesquelles on cherche à faire l'interpolation. Voici comment cela se passe :

  • A un instant t1, on a fait un calcul MCNP d'où résultent des < $ \sigma$$ \varphi$ > 1.
  • REM lance un calcul IME à < $ \sigma$$ \varphi$ > 1 constant, (seul le paramètre d est non nul). Il demande l'intégration jusqu'à un intant t2 (pas MCNP) pour lequel IME trouve une composition $ \overline{{N_{2}}}$ (la mémorisation des résultats intermédiaires - ajustements, synchronisations s'effectue pendant cet intervalle).
  • REM lance un calcul MCNP pour l'instant t2 avec la composition $ \overline{{N_{2}}}$, ce qui donne les < $ \sigma$$ \varphi$ > 2 recherchés.
  • REM calcule les paramètres du spline d'interpolation entre < $ \sigma$$ \varphi$ > 1 et < $ \sigma$$ \varphi$ > 2 et relance le calcul IME entre les instants t1 et t2 à < $ \sigma$$ \varphi$ > variable (il passe les paramètres a, b, c, d du fit du spline dans la matrice d'intégration). Il en résulte une nouvelle composition pour l'instant t2, soit $ \overline{{N_{2}}}{^\prime}$.
  • REM vérifie si l'écart quadratique absolu entre $ \overline{{N_{2}}}$ et $ \overline{{N_{2}}}{^\prime}$ [ $ \chi^{{2}}_{}$ = ($ \overline{{N_{2}}}$ - $ \overline{{N_{2}}}{^\prime}$)2] est inférieur à une consigne fournie par le programme principal.
  • Si oui, la composition $ \overline{{N_{2}}}{^\prime}$ est conservée. Le calcul IME reprend à partir de t2 jusqu'à un temps t3 tel que ( t3 - t2) = (t2 - t1), avec la composition $ \overline{{N_{2}}}{^\prime}$ et < $ \sigma$$ \varphi$ > 2 constant, afin de trouver < $ \sigma$$ \varphi$ > 3, interpoler, etc.
  • Dans le cas où le $ \chi^{{2}}_{}$ est très bon (inférieur à une deuxième consigne), le temps total d'intégration fixé à IME (le pas MCNP) peut être allongé. REM spécifie alors t3 tel que ( t3 - t2) > (t2 - t1).
  • Si, par contre, l'écart entre $ \overline{{N_{2}}}$ et $ \overline{{N_{2}}}{^\prime}$ est supérieur à la consigne, un nouveau calcul MCNP est lancé pour l'instant t2, avec la composition $ \overline{{N_{2}}}{^\prime}$, ce qui donne un < $ \sigma$$ \varphi$ > 2'.
  • REM calcule les paramètres du spline d'interpolation entre < $ \sigma$$ \varphi$ > 1 et < $ \sigma$$ \varphi$ > 2' et relance le calcul IME entre les instants t1 et t2 à < $ \sigma$$ \varphi$ > variable . Il en résulte une nouvelle composition pour l'instant t2, soit $ \overline{{N_{2}}}{^\prime}{^\prime}$.
  • Si l'écart quadratique entre $ \overline{{N_{2}}}{^\prime}$ et $ \overline{{N_{2}}}{^\prime}{^\prime}$ est inférieur à la consigne, on garde $ \overline{{N_{2}}}{^\prime}{^\prime}$. Le calcul IME reprend à partir de t2 jusqu'à un temps t3 tel que ( t3 - t2) = (t2 - t1), avec la composition $ \overline{{N_{2}}}{^\prime}{^\prime}$ et < $ \sigma$$ \varphi$ > 2' constant, afin de trouver < $ \sigma$$ \varphi$ > 3, interpoler, etc.
  • Si l'écart quadratique entre $ \overline{{N_{2}}}{^\prime}$ et $ \overline{{N_{2}}}{^\prime}{^\prime}$ est supérieur à la consigne, on réduit t2 : on diminue l'intervalle t2 - t1 et on recommence l'intégration IME à partir de t1, en procédant comme ci-dessus pour les vérifications de précision, et en allongeant le pas MCNP dès que possible, c'est-à-dire quand on trouve un très bon écart quadratique entre $ \overline{{N_{k}}}$ et $ \overline{{N_{k}}}{^\prime}$.

Avec ce système de contrôle de précision et d'allongement du pas MCNP dès que c'est possible, la fréquence des calculs MCNP est adaptée au comportement du réacteur. Si le réacteur est dans une phase transitoire, les pas MCNP sont courts, comme le sont les pas d'intégration d'IME, comme le sont aussi les temps d'ajustement. Lorsque le réacteur tend vers l'équilibre, le temps entre deux calculs MCNP augmente, les durées d'intégration IME sont de plus en plus longues. Cependant, grâce aux ajustements et synchronisations imposées à IME, on dispose de mesures intermédiaires qui rendent compte de l'évolution des compositions, de la puissance, de la réactivité, on dispose de la possibilité de modifier des paramètres du réacteur.

La façon dont s'intercalent les différentes étapes de calcul IME et MCNP est illustrée dans la figure suivante.

integration_REM.jpg
Quelques étapes du calcul d'évolution d'un réacteur montrant que les pas d'ajustement, les pas MCNP et les durées d'intégration RK4 changent au cours du temps, et pas forcément de façon régulière.


Simplification de l'intégration

Certains noyaux du vecteur composition ont une durée de vie suffisamment courte pour qu'on puisse considérer qu'ils sont en équilibre avec leur père quand le réacteur a fonctionné assez longtemps pour que leur nombre ait eu le temps de se stabiliser avec la population du père ; on peut accélérer le calcul en les mettant en équilibre avec le père, c'est-à-dire qu'ils restent dans le vecteur $ \overrightarrow{N} $ mais leur population, au lieu d'être calculée en fonction des éléments de la matrice d'intégration, est fixée en fonction de la population du père.

Prenons comme exemple le 135Xe dans un réacteur à sels fondus. C'est un gaz, il sort rapidement du coeur du réacteur grâce à un bullage d'hélium. Sa période dans le coeur est de 30 secondes. Cependant, il a une section efficace de capture grande, si bien que son taux de capture, n < $ \sigma$$ \varphi$ > ne peut être négligé : les captures sur le Xenon font perdre des neutrons proportionnellement à la quantité de Xenon présent, ne pas en tenir compte changerait la réactivité du réacteur. On est donc obligé de suivre le Xenon. Mais, si on le suivait comme les autres noyaux, les bruits de calcul introduiraient des fluctuations du vecteur $ \overrightarrow{N} $ telles qu'on ne parviendrait jamais à allonger le pas d'intégration dans IME. D'où ce compromis, qui consiste à estimer la population du noyau en fonction de la population de son père, quand c'est possible, d'autant que la population peut être mieux estimée par l'hypothèse de l'équilibre avec le père que par le calcul. Bien entendu, le Xenon n'est pas le seul noyau qui peut être traité ainsi.

C'est IME qui assure la mise à l'équilibre d'un noyau avec son père, sous le contrôle de REM. Celui-ci donne à IME un temps, le seuil de demie vie effective (calculée en combinant le taux de réaction et la décroissance radioactive du noyau) en dessous duquel IME met le noyau en équilibre avec son père. Cependant, il faut être sûr que le réacteur a fonctionné assez longtemps pour que cette mise à l'équilibre soit valable. C'est pourquoi REM donne à IME un temps minimal de fonctionnement, à partir duquel la mise à l'équilibre est autorisée.


La boucle de calcul

En résumé, la boucle de calcul est une suite de calculs MCNP, suivis de calculs , le tout étant piloté par la précision, comme illustré dans la figure. Dans REM, c'est le « pas MCNP » (qui est aussi la durée totale d'intégration IME) qui est ajusté en fonction de la variabilité des vecteurs $ \overrightarrow{N} $.

REM reste dans sa boucle de calcul jusqu'à ce que le temps total de simulation du réacteur, fixé par le programme principal, soit atteint et là, il lui rend la main. Notons que le programme principal a pu suivre le processus grâce à l'activation de ses procédures en cours de calcul, lors des ajustements et/ou synchronisations).


Gestion des compositions et des cellules volumiques

Voyons comment sont gérés les éléments de volume tels qu'ils sont vus par REM, IME, CGM.

Conceptuellement, un réacteur se découpe en zones dont la fonction, donc le type de matériaux qu'elles contiennent, diffèrent. On est amené, par exemple, à caractériser le combustible, le modérateur, le caloporteur, les éléments de structure, les réflecteurs. Ces matériaux, que nous appellerons des « composites » sont à des températures différentes suivant leur emplacement dans le réacteur. Pour un composite donné, par exemple un réflecteur, la composition initiale est la même, quelle que soit la température de fonctionnement. Par contre, elle évolue de façon différente selon la zone de température.

Par ailleurs, les éléments de volume, les cellules au sens CGM, que nous appellerons les « cellules CGM », qui décrivent le découpage géométrique du réacteur demandé par le programme principal, peuvent fractionner une zone de température. La correspondance entre les « cellules CGM » et les zones de température doit être connue de REM, elle est fournie par le programme principal.

Ainsi, REM reçoit du programme principal une liste de matériaux - les composites - avec leur composition initiale, les zones de température (et leur température), la liste des cellules CGM correspondantes.

REM construit les arbres d'évolution pour chaque composite par zone de température, et non pas pour chaque cellule CGM, ce qui serait inutile : on obtiendrait une duplication d'arbres d'évolution dans le cas où plusieurs cellules CGM recouvrent une même zone de température.

Les calculs MCNP sont toujours faits individuellement sur les cellules CGM. Par contre, les calculs IME peuvent être faits sur un groupe de cellules qui recouvre tout ou partie d'une zone de température pour un composite donné. Un tel groupement de cellules CGM suppose qu'elles ont la même composition à tout instant (calcul IME), bien qu'elles ne soient pas forcément soumises exactement au même flux (calcul MCNP).

Une fois que les arbres d'évolution sont construits, REM peut déterminer quels sont les tallys qu'il voudra passer à CGM pour qu'il les demande à MCNP (ils figurent dans le fichier de commande MCNP). En ce qui concerne les sections efficaces, c'est simple, c'est les < $ \sigma$$ \varphi$ > de toutes les réactions possibles, pour chacun des noyaux susceptible d'interagir (pour lesquels ACE a trouvé des sections efficaces) qui figurent dans l'arbre d'évolution. En effet, on démontre que MCNP peut calculer les sections efficaces moyennes d'interactions sur des noyaux qui ne sont qu'à l'état de trace (ou pas encore apparus dans la composition fournie à MCNP) du moment que les fichiers de sections efficaces lui sont fournis. Il suffit que MCNP ait évalué le flux (qui n'est, a priori, pas influencé par des noyaux qui ne sont dans le milieu qu'à l'état de trace). Ainsi, il est possible de suivre l'évolution des noyaux qui ne figurent pas dans les compositions initiales au fur et à mesure qu'ils apparaissent : leurs < $ \sigma$$ \varphi$ > étant évalués, ils participent à la matrice d'évolution passée à IME.


Décroissances fictives

Lors de la simulation d'un réacteur, on peut être amené à simuler l'extraction de certains noyaux : dans un RSF, par exemple, on utilise un bullage d'helium pour extraire des produits de fission gazeux. Pour cela, nous avons introduit la notion de décroissance radioactive fictive. Soit un noyau i, bien réél. Il a ses sections efficaces d'interaction (fournies par ACE) et/ou ses décroissances radioactives (fournies par DPH). Lui ajouter une décroissance fictive, c'est lui ajouter une décroissance radioactive vers un noyau, disons k, « fictif », c'est-à-dire, en l'occurence, un noyau qui n'interagira pas avec les neutrons. En effet, un noyau extrait est, par définition hors flux, il ne peut plus interagir avec les neutrons. Par contre, on veut suivre son évolution et, s'il est radioactif, on veut pouvoir suivre ses descendants qui, eux aussi, doivent être hors flux.

L'artifice que nous avons utilisé est d'attribuer un numéro d'isomère négatif aux noyaux « fictifs », tel que leur numéro d'isomère d'origine puisse être retrouvé et, par conséquent, leur éventuelle décroissance radioactive suivie (qui se fera vers un noyau dont l'isomère est négatif, puisqu'il faut rester hors flux). ACE n'est pas interrogé pour les noyaux dont l'isomère est négatif, ils n'ont aucune section efficace d'interaction. Mais on veut suivre leur évolution avec IME : sachant quels noyaux il faut « extraire » dans chaque composite, et avec quelle constante de temps chacun « décroît » vers le noyau fictif associé, REM ajoute ces noyaux (et leurs descendants - fictifs également) aux arbres d'évolution. Ils sont ensuite traités de façon standard.

Ce dont on dispose à la fin de REM

On a déjà vu que les compositions trouvées à chaque ajustement et synchronisation lors du calcul d'évolution sont mémorisées, pour toutes les cellules que connaît IME, qui peuvent correspondre à des groupes de cellules CGM. En outre, on mémorise tous les tallys demandés à MCNP. Afin d'en faciliter l'exploitation, ces informations sont organisées de manière à regrouper en un ensemble

  • les résultats de l'évolution durant tout un pas MCNP
  • les résultats du calcul MCNP de fin de ce pas

Nous appellerons un tel ensemble un inventaire.

Pour chaque inventaire, pour chaque cellule, on a donc mémorisé :

  • les compositions de chaque instant caractéristique (ajustement, synchronisation)
  • les tallys déterminés par MCNP
  • la réactivité telle que déterminée par MCNP (qui la calcule différemment de IME) et, pour chaque instant, des données globales du réacteur, soit :
  • le nombre de neutrons créés
  • le nombre de neutrons absorbés
  • la réactivité évaluée par IME
  • le nombre de neutrons source (le paramètre $ \Phi$ sur lequel joue IME pour ajuster la puissance)

spline : courbe du 3 eme degré couramment utilisée pour l'interpolation entre des points expérimentaux. Elle assure la continuité des dérivées aux points de jonction entre deux splines (points expérimentaux).

La méthode Monte Carlo consiste à suivre l'histoire de chaque neutron dans le réacteur, de sa « naissance » (source externe, neutrons de fission, ...) à sa « mort » (capture par un noyau ou fuite hors du réacteur). L'histoire de chaque neutron dépend des interactions qu'il a avec la matière. Le type de noyau sur lequel a lieu l'interaction est tiré au hasard (comme dans les jeux de hasard du Casino de Monte Carlo). Il dépend d'une part de la section efficace de chaque noyau, d'autre part du nombre relatif de noyaux du type considéré par rapport au nombre total de noyaux. Ensuite, le type d'interaction est choisi, lui aussi au hasard, en fonction des sections efficaces relatives de chaque type d'interaction possible à l'énergie du neutron. Ainsi, on obtient un comportement moyen du réacteur à un instant donné (flux, coefficient de criticité, puissance, ...) en suivant de nombreuses histoires individuelles de neutron (de 100 000 à 1 000 000 de neutrons).

Pour lancer un calcul MCNP, il faut fournir :

Une description de la géométrie du réacteur aussi détaillée que possible.

La composition précise de chaque matériau (corps présents, densité, ...).

La température dans les diverses parties : la dilatation des matériaux en dépend, et celle-ci modifie la distance parcourue par les neutrons entre les interactions (le libre parcours moyen). Une augmentation de température induit également des variations sur les sections efficaces, comme l'élargissement apparent des résonances (effet Doppler).

Les chemins d'accès aux fichiers où se trouvent les sections efficaces d'interaction en fonction de l'énergie des neutrons.

Les tallys, c'est-à-dire les quantités qu'il faut calculer (flux, réactivité, sections efficaces moyennes des réactions, ...).