Méthodes de calcul neutronique

de cœur
par Giovanni B. BRUNA
Docteur en Physique Nucléaire, Université de Gênes
Chargé de missions au Département Performances de Cœur de la société Framatome
et Bernard GUESDON
Ingénieur de l’École Nationale Supérieure de Mécanique et du Génie Atomique
Adjoint technique au Chef de Division Procédé de la société Framatome

1. Démarche de conception neutronique de cœur.............................. B 3 070 - 2

2. Techniques de calcul............................................................................... — 4
2.1 Équation de Boltzmann et approximation de la diffusion........................ — 4
2.2 Séparation des étapes : calculs de cellule et calculs de cœur ................. — 5
2.3 Génération des sections efficaces pour les calculs de cœur ................... — 6
2.4 Approximations numériques...................................................................... — 7
2.5 Corrections ponctuelles des opérateurs .................................................... — 8
3. Calculs du transport des neutrons...................................................... — 8
3.1 Méthodes déterministes (analytiques) ...................................................... — 8
3.2 Méthodes probabilistes (Monte-Carlo) ...................................................... — 10
4. Méthodes de résolution de l’équation de la diffusion .................. — 12
4.1 Méthodes aux différences et aux éléments finis ...................................... — 12
4.2 Discrétisation et effet de maillage.............................................................. — 12
4.3 Processus de calcul ..................................................................................... — 13
4.4 Domaines d’application .............................................................................. — 14
5. Méthodes nodales avancées ................................................................. — 15
5.1 Principe de l’équivalence nodale................................................................ — 15
5.2 Méthode d’expansion nodale ..................................................................... — 16
5.3 Prise en compte de l’épuisement ............................................................... — 17
5.4 Reconstruction fine...................................................................................... — 17
5.5 Domaines d’application .............................................................................. — 18
6. Méthodes de perturbation en physique des réacteurs.................. — 19
6.1 Méthode CPT (Classical Perturbation Theory ).......................................... — 19
6.2 Méthode GPT (Generalized Perturbation Theory ).................................... — 20
6.3 Développements récents ............................................................................ — 20
7. Schémas de calcul ................................................................................... — 21
7.1 Structure et qualification des chaînes de calcul........................................ — 21
7.2 Chaîne de calcul de conception des REP ................................................... — 22
7.3 Formulaire de projet des RNR .................................................................... — 23
Annexe — Équation générale du transport ............................................... — 24
2 - 1996

Pour en savoir plus........................................................................................... Doc. B 3 070

es calculs neutroniques de cœur conduisent à la détermination de la géo-

L
B 3 070

métrie, de la composition, des moyens de contrôle et des limites d’exploi-

tation du cœur en accord avec les performances attendues, les règles de
sûreté, les options technologiques et les modes de gestion du combustible.

Toute reproduction sans autorisation du Centre français d’exploitation du droit de copie est strictement interdite.
© Techniques de l’Ingénieur, traité Génie nucléaire B 3 070 − 1
MÉTHODES DE CALCUL NEUTRONIQUE DE CŒUR _____________________________________________________________________________________________

Ce chapitre présente les méthodes actuelles de calcul neutronique appliquées

par l’industrie nucléaire. Ces méthodes se caractérisent par l’utilisation d’outils
de calcul modernes et par une validation étendue dans un cadre formel.
Le paragraphe 1 décrit l’objectif des calculs neutroniques et présente, en les
justifiant physiquement, les hypothèses et les approximations qui interviennent
dans la démarche qui conduit à l’équation du bilan (le développement est
présenté en Annexe). Sont, en outre, décrites et définies mathématiquement la
réactivité et la distribution de puissance, grandeurs qui, dans des domaines
différents, jouent un rôle primordial dans l’exploitation des réacteurs nucléaires.
Dans le paragraphe 2 sont présentés les éléments de base de la technique de
calcul neutronique de cœur. On y aborde, notamment, les problèmes essentiels
de la discrétisation de l’équation du bilan neutronique du système en espace
et en énergie, et de la prise en compte et du traitement des hétérogénéités éner-
gétiques et spatiales en milieu multiplicateur. On explique également comment,
dans les chaînes modernes, ces deux aspects sont traités séparément, dans les
codes de calcul dits de cellule et de cœur. On y traite, enfin, de la prise en
compte de la variable temporelle en fonctionnement nominal, non accidentel.
Les paragraphes 3, 4 et 5 présentent un panorama des méthodes de calcul
utilisées dans les différents domaines de la conception des réacteurs à neu-
trons rapides et thermiques. Les aspects plus strictement numériques y
côtoient les considérations d’intérêt général relatives au domaine d’utilisation
des différentes méthodes et aux options adoptées industriellement.
Le paragraphe 6 présente un panorama complet des utilisations industrielles
des méthodologies dérivées de la théorie des perturbations, dans ses formes
classique et généralisée.
Le paragraphe 7 traite le problème crucial de la qualification et de la définition
du domaine de validité des chaînes de calcul. Les deux exemples présentés font
directement appel à l’expérience industrielle acquise en France dans les réac-
teurs à eau préssurisée (REP) et dans les réacteurs à neutrons rapides (RNR).
Avant de commencer ce chapitre, il est conseillé au lecteur de se reporter aux
chapitres suivants dans ce traité :
— le chapitre Physique du cœur : neutronique et thermohydraulique
[B 3 090], qui développe le dimensionnement et la détermination des caracté-
ristiques physiques des réacteurs intégrant l’ensemble des critères de
conception thermiques, neutroniques, technologiques et économiques ;
— les chapitres Thermohydraulique des réacteurs [B 3 050] et Théorie des
réacteurs nucléaires [B 3 025], qui développent la thermohydraulique et la neu-
tronique.
Nota : ce chapitre est une œuvre collective de compilation à laquelle ont contribué plusieurs
spécialistes, chacun dans son domaine de compétence. Nous remercions pour leur collabora-
tion Messieurs Aldo Dall’Osso, Richard Aigle, Michel Doucet, Patrick Girieud de Framatome et
Giovani Gastaldo de Novatome.

1. Démarche de conception — de l’imprécision résultant des approximations numériques et

des simplifications du calcul ;
neutronique de cœur — du caractère prévisionnel des études.
Les plages d’incertitude imputables aux données de base et aux
méthodes de calcul ressortissent de la qualification ; elles sont
Au stade des études de conception, la définition des limites imposées au concepteur qui doit, en revanche, prendre en compte
d’exploitation s’accompagne de l’évaluation et de la prise en la variabilité des autres paramètres (tolérances de fabrication,
compte d’incertitudes et, éventuellement, de marges assurant le réponses des instruments de mesure, modes de gestion, etc.) et
conservatisme des résultats, en raison : aléas.
— de la connaissance imparfaite que l’on peut avoir de certaines L’état du cœur d’un réacteur nucléaire est régi par deux paramè-
données de base (sections efficaces, propriétés des matériaux, tres principaux (cf. chapitre Théorie des réacteurs nucléaires
etc.) ; [B 3 025]), à savoir :

Toute reproduction sans autorisation du Centre français d’exploitation du droit de copie est strictement interdite.
B 3 070 − 2 © Techniques de l’Ingénieur, traité Génie nucléaire
____________________________________________________________________________________________ MÉTHODES DE CALCUL NEUTRONIQUE DE CŒUR

— la réactivité, qui mesure en transitoire l’écartement de la * L’unicité de la solution δΦ = 0 est garantie par l’alternative de Fredholm (cf. chapitre
Équations aux dérivées partielles [A 650] dans le traité Sciences fondamentales) car le
population neutronique par rapport aux conditions d’équilibre et, terme de source défini par l’équation (2) n’est pas orthogonal à la solution de l’adjointe
en conditions stationnaires, le défaut de bilan neutronique. S’agis- de l’équation (1).
sant d’une grandeur intégrale, ce paramètre est, à tout instant et On transforme ainsi une équation, inhomogène par nature, en
pour une configuration de cœur donnée, unique et indépendant de une équation homogène, d’où un gain de simplicité et, surtout, de
la position ; rapidité de calcul**.
— la puissance, qui traduit en termes macroscopiques la distri-
bution spatiale des événements de fission dans le cœur. Ce para- ** En réalité, la forme générique du terme  est :  = ρ P 0 Φ 0 + δ  , où la composante
mètre dépend des caractéristiques des milieux et de la distribution δ est définie dans le sous-espace orthogonal à Φ 0 . Dans le cas général, l’équation (4)
en espace et en énergie des neutrons à l’intérieur du système est remplacée par le système :
multiplicateur. Il peut être, de ce fait, extrêmement sensible aux
 A0 Φ0 + P0 ( 1 – ρ ) Φ0 = 0
conditions locales.  (7)
 A 0 δΦ + P 0 δΦ = δ 
Pour le concepteur, l’objectif des calculs neutroniques de cœur
est d’évaluer ces deux paramètres avec la précision maximale δΦ, composante de la solution orthogonale au fondamental, a alors valeur non nulle et
compatible avec la taille et la complexité des problèmes à traiter, le flux conserve toutes ses composantes, y compris celles non calées sur le mode fonda-
mental.
l’utilisation de méthodes numériques et les limites imposées à leur
qualification par les incertitudes expérimentales. Les grandeurs Cette procédure mathématique mérite quelques commentaires :
ainsi calculées, assorties d’incertitudes et de marges, sont — en imposant la proportionnalité entre le terme inhomogène et
comparées aux valeurs limites de conception établies à partir des le flux neutronique [relation (2)], on vide arbitrairement celui-ci de
critères de sûreté, pour vérifier le respect de ces derniers. tous les modes autres que le fondamental : en d’autres termes, on
La réactivité et la distribution de puissance, qui revêtent une impor- projette le système sur le mode fondamental en engendrant arti-
tance égale dans les études de conception, ne sont toutefois pas des ficiellement, entre les valeurs réelle et calculée du flux, des
grandeurs homogènes et ne relèvent pas du même domaine de défi- différences qui pourraient être notables dans quelques cas parti-
nition. culiers (réacteur sous-critique à l’arrêt, transitoires, etc.) ;
— en revanche, la réactivité est, de par sa définition, tout à fait
Alors que la puissance est, à tout instant, l’image de la distribu- insensible aux modes supérieurs : liée au mode fondamental, sa
tion instantanée des fissions dans le système multiplicateur (elle valeur est identique si on la calcule dans l’état réel du système ou
représente donc une réalité objective), la réactivité est une grandeur dans sa projection sur le mode fondamental. On dispose ainsi d’un
qui, en conditions stationnaires, ne possède pas de véritable signi- paramètre absolu qui ne dépend que des propriétés intrinsèques
fication physique. Objet mathématique, elle traduit l’exigence, tou- du système multiplicateur et reste insensible à son état.
jours présente dans les études, de décrire, par des équations
statiques et homogènes, des situations dans lesquelles la popula- L’équation (4) peut aussi s’écrire sous d’autres formes.
tion neutronique ne serait pas en équilibre. Dans ce cas, au lieu En définissant :
d’introduire explicitement dans le bilan, qui lie le nombre total de 1 1
disparitions de neutrons par absorption et/ou fuite A0 Φs au nombre k eff = -------------- et λ 0 = ---------- (8)
1–ρ k eff
de production par fission P0 Φs , un terme inhomogène relevant de
ce déséquilibre : où le paramètre k eff prend le nom de facteur de multiplication
A 0 Φ s + P0 Φ s =  (1) effective, on obtient les formes canoniques :
avec A0 et P0 opérateurs d’absorption plus fuites (disparition) et P0 
 A + --------
- Φ = 0 (9)
de production du système multiplicateur,  0 k eff 0
Φs flux de neutrons,
on fait implicitement l’approximation : ou encore :
A 0 Φ 0 + λ0 P 0 Φ 0 = 0 (10)
 = ρ P0 Φ 0 (2)
L’équation (10) admet une adjointe dans laquelle figurent les
avec : t t
opérateurs A 0 et P 0 , transposés de A 0 , P0 :
Φs = Φ0 + δΦ (3)
t
où ρ est une constante de proportionnalité qui prend le nom A t0 Φ*0 + λ 0 P 0 Φ*0 = 0 (11)
de réactivité,
Φs est la solution de l’équation inhomogène (1), dont la solution Φ*0
, est appelée flux adjoint ou importance neu-
tronique. Cette définition n’est strictement rigoureuse que lorsque
δΦ est la différence entre Φs et Φ0 ’ Φ0 (mode fondamental) t t
les opérateurs A 0 et P 0 coïncident avec les inverses de A0 et P0 .
étant la solution de l’équation aux valeurs propres asso-
ciée à l’équation (1) : En multipliant les deux membres de l’équation (2) par Φ *0 et en
réarrangeant les termes, on obtient la définition générique de la
réactivité :
A 0 + ( 1 – ρ )P 0 Φ0 = 0 (4)
〈 Φ* 0 , 〉
ρ = ----------------------------------- (12)
Cette approximation équivaut à prendre en compte la seule 〈 Φ* 0 ,P 0 Φ 0〉
composante de Φs parallèle au mode fondamental, car la fonction
δΦ doit satisfaire strictement l’équation : où le symbole < > indique, classiquement, une intégration sur
l’espace (ici la totalité du système) et sur l’énergie.
(A0 + P0) (Φ0 + δΦ) = ρ P0 Φ0 (5) Nota : pour plus d’information sur la définition de la réactivité, se reporter à [2] [17]
et au paragraphe 6.
qui, tenant compte de l’équation (4), implique la relation :
(A0 + P0) δΦ = 0 (6)
laquelle, en conformité avec l’alternative de Fredholm* [1], admet
comme seule et unique solution : δΦ = 0.

Toute reproduction sans autorisation du Centre français d’exploitation du droit de copie est strictement interdite.
© Techniques de l’Ingénieur, traité Génie nucléaire B 3 070 − 3
MÉTHODES DE CALCUL NEUTRONIQUE DE CŒUR _____________________________________________________________________________________________

2. Techniques de calcul et :

νΣ f, 1 0 . 0
0 νΣ f, 2 . 0
2.1 Équation de Boltzmann [ νΣ f, gg′ ] =
0 0 . 0
et approximation de la diffusion 0 0 . νΣ f, G

En conditions stationnaires, le bilan neutronique d’un système où [χg, g ’] este la matrice de distribution en énergie des neutrons
nucléaire multiplicateur est régi par une équation, en général issus de fission et [ν ∑f, gg’ ] l’opérateur de fission multigroupe du
inhomogène. Moyennant des hypothèses simplificatrices, cette système.
équation peut, dans un bon nombre de cas pratiques, être écrite Le couplage en énergie entre les neutrons qui engendrent les
sous la forme générique d’une équation aux valeurs propres fissions et ceux qui en sont issus est négligé, en s’appuyant sur
appelée équation de Boltzmann (10). une analyse statistique. Cela permet de simplifier le terme [χg, g’]
La résolution de cette équation par voie numérique comporte sa en le transformant en un vecteur de la forme :
discrétisation en espace, en énergie et en probabilité d’interaction
et/ou de direction des neutrons, selon l’approximation choisie pour χ1
l’opérateur de transport (fuites, cf. Annexe). Or, le système d’équa- χ2
tions obtenu après discrétisation reste en général trop complexe [ χg ] =
pour les performances de calcul des ordinateurs utilisés. Il est alors .
résolu en faisant appel à d’autres simplifications. C’est ainsi que, χG
dans les études, on adopte l’approximation de la diffusion.
Dans la théorie de la diffusion, l’équation d’état du système peut Les opérateurs A0 et P0 , qui ne dépendent pas explicitement du
être obtenue à partir de celle du transport, soit sur la base de temps, sont homogènes dans des sous-domaines en espace et en
considérations physiques [3], soit mathématiquement par un déve- énergie, et sont variables en dimension selon la discrétisation
loppement asymptotique (par exemple [18]). choisie.
L’écriture de l’équation de la diffusion comporte plusieurs Écrite sous forme explicite pour chaque maille i et pour chaque
approximations et simplifications par rapport à la forme intégrale groupe d’énergie g, l’équation (13) devient :
d’origine de l’équation du bilan en transport :

L
J dS œ +  Σ a, i + ∑ Σ i
g → g′
— le flux est supposé isotrope en tout point de l’espace et à
∑
g g g
+ Φoi V i (14)
toute énergie, moyennant toutefois une correction d’anisotropie  
œ=1 Sœ g′
introduite au niveau du bilan comme terme correctif global du
–∑  Σ i + λ 0 χ ν Σ f, i Φoi V i = 0
nombre d’interactions entre les neutrons et la matière ; g′ → g g g′ g′
— la géométrie est simplifiée par homogénéisation, notamment  
g′
en ce qui concerne la description des crayons et parfois des
assemblages ; avec i et g indiquant respectivement une maille et un groupe
— le courant est approximé en tous points au premier ordre en d’énergie, L le nombre de mailles adjacentes à i, dS œ la surface de
g
imposant la loi de Fick [3] : J0 = – D grad Φ0 , avec D coefficient de contact entre i et une des mailles génériques adjacentes œ, Φ 0i le
diffusion rendant compte des propriétés des milieux. flux moyen du groupe g et Vi le volume de la maille i.
L’équation générique de la diffusion s’écrit, après condensation Pour les calcules des REP, on retient habituellement deux groupes
en énergie, homogénéisation en espace et discrétisation des opé- d’énergie : le premier, dit rapide, englobe tous les neutrons d’éner-
rateurs en espace et énergie : gie supérieure à 0,625 eV, le second, dit thermique, regroupe les
neutrons d’énergie inférieure à 0,625 eV (1,82 eV chez certains
 A 0AΦ0 + λ0 P
0 0 Φ0 =
. 0 . 0  (13) concepteurs). Cette répartition est justifiée par la dépendance de
 1 
 l’énergie des sections efficaces de fission de l’U235, utilisé généra-
avec λ0 = 1/keff –Σ2 → 1 A2 0 . 0 
  lement comme combustible dans les REP. L’augmentation exponen-
où : A0 =  –Σ3 → 1 –Σ3 → 2 Ag . .  tielle des sections efficaces de l’U235 dans le domaine thermique
  explique que, dans ce type de réacteurs, les neutrons appartenant
 . . . . . 
 . –Σ  au deuxième groupe, en dépit de leur nombre limité, engendrent
 –Σ G→1 –Σ G→2 G→G–1 AG  plus de 70 % des fissions [3] [4].
Pour les calculs des RNR, on adopte en général un découpage à
et : Ag = – div [D (r ) grad] + Σag (r ) 1  g  G
six groupes d’énergie centré sur les domaines rapide et épithermi-
que du spectre (on entend par domaine rapide la région d’énergie
P0 = [χg, g’] [ν Σf, gg’] 1  g, g′  G
supérieure à 500 keV, par domaine épithermique celle comprise
G étant le nombre de groupes énergétiques, Σ la section efficace, entre 500 keV et quelques eV où se situent les grandes résonances
r le vecteur position, avec de façon générique : des sections efficaces neutroniques). Un découpage plus fin, en 25
groupes, est utilisé si besoin en fonction des problèmes traités [3]
[4].
χ 1, 1 χ 2, 1 . χ G, 1
En ce qui concerne l’espace, les calculs sont généralement effec-
χ 1, 2 χ 2, 2 . χ G, 2
[ χ g, g′ ] = tués au niveau du crayon, car c’est ponctuellement qu’il faut véri-
. . . . fier le respect des limites de conception et des critères de sûreté.
χ 1, G χ 2, G . χ G, G

Toute reproduction sans autorisation du Centre français d’exploitation du droit de copie est strictement interdite.
B 3 070 − 4 © Techniques de l’Ingénieur, traité Génie nucléaire
____________________________________________________________________________________________ MÉTHODES DE CALCUL NEUTRONIQUE DE CŒUR

Selon la géométrie du problème, l’équation (13) est assortie de Dans les méthodes de conception, on partage le calcul en deux
conditions aux frontières : flux ou courant nul aux interfaces. Elles parties, ce qui permet de traiter indépendamment les deux
sont de la forme : niveaux d’hétérogénéité :
∂ Φ0 — d’une part, l’interaction des neutrons avec la matière, avec ses
α Φ 0 + β ---------- = 0 (15) spécificités énergétiques et spatiales qui relèvent d’effets locaux liés
∂n
à la géométrie des milieux, à la dépendance énergétique des sec-
où n indique une direction générique et α, β des constantes de pro- tions efficaces et à la distribution des neutrons en énergie (spectre
portionnalité. neutronique). Ce premier aspect comporte une modélisation très
Le système a une solution unique. Cela peut être démontré en précise et détaillée des milieux physiques à un nombre élevé de
–1 groupes d’énergie, mais limitée à une région restreinte en espace
définissant l’opérateur intégral Q = A 0 P 0 et en appliquant le théo- (de l’ordre du libre parcours moyen des neutrons). Il est traité par
rème de Perron-Frobenius [1] [5]. Le flux, solution de l’équation, est les codes de cellule, appelés aussi codes de spectre, dont le rôle
alors déterminé à un facteur multiplicatif près. principal est celui de générer des sections efficaces pour les calculs
de cœur (§ 2.3) ;
En général, on choisit comme condition de normalisation la puis-
— d’autre part, le bilan neutronique du système, paramètre plus
sance nominale du réacteur, exprimée par l’intégrale 〈 κΣ f , Φ0 〉, où
global et macroscopique qui traduit, en premier abord, le couplage
le paramètre κ indique l’énergie thermique dégagée en moyenne
entre régions du cœur, couplage qui est entretenu majoritairement
dans le système pour chaque fission (*), ou bien on impose une
par les neutrons rapides peu sensibles aux effets locaux. Cet
source de fission unitaire.
aspect nécessite une description du système dans son ensemble
(*) L’énergie directement dégagée dans le cœur provient du ralentissement des frag-
ments de fission, de l’atténuation des γ de fission, du ralentissement des neutrons issus
qui peut être très simplifiée en ce qui concerne le nombre de grou-
de fission, de l’atténuation des γ de capture, de l’excitation des noyaux par capture neu- pes d’énergie et la prise en compte de l’hétérogénéité locale des
tronique. Le terme κ tient compte, de façon globale et forfaitaire, de tous ces apports, milieux. Cette description est du domaine des codes de cœur.
ramenés aux sources de fission. Dans le cas d’un combustible standard, la contribution
des captures s’élève à quelques pourcent du total, mais cette valeur augmente cependant Les codes de cellule présentent des particularités qui tiennent de
sensiblement en présence de forts absorbants, poisons consommables et/ou grappes de la sophistication des modèles physiques et du type de données à
contrôle. fournir aux codes en aval (on verra par la suite comment ces don-
nées dépendent étroitement du mode de résolution adopté).
Toutefois, tous ont des caractéristiques communes, à savoir :
2.2 Séparation des étapes : calculs — ils utilisent une quantité extrêmement importante de données
de base (sections efficaces microscopiques de tous les isotopes
de cellule et calculs de cœur intervenant dans la composition des milieux, concentrations isoto-
piques, géométrie détaillée, etc.) ;
— ils font appel à un ensemble d’approximations et de modéli-
Afin de rechercher la solution de l’équation (13) assortie de ses
sations physiques qui s’appuient sur un processus de calcul numé-
conditions aux frontières (15), il est indispensable de connaître
rique sophistiqué, dans lequel elles sont étroitement imbriquées ;
explicitement la valeur de toutes les composantes spatiales et
— ils gèrent des masses énormes d’informations non seulement
énergétiques des opérateurs matriciels A0 et P0 , dans chaque élé-
pendant l’exécution du calcul proprement dit, mais aussi en sortie
ment de volume de l’espace (hyper-espace) dans lequel le pro-
quand il s’agit de mettre en forme les données à transmettre aux
blème est défini et borné. Dans l’approximation de la diffusion, la
codes de cœur.
dimension de cet espace est égale au nombre de dimensions spa-
tiales considérées plus le nombre de groupes. Les codes de cœur diffèrent assez considérablement entre eux
en ce qui concerne les options retenues pour discrétiser l’équation
Or en raison des spécificités physiques du processus d’interac-
du bilan neutronique et les méthodologies adoptées dans sa réso-
tion des neutrons avec la matière, de sa dépendance en énergie et
lution, mais ils restent très proches sur le plan de la modélisation
de la géométrie des milieux, le calcul de ces composantes est, en
physique et des résultats fournis. Seule la prise en compte de
général, très complexe. Pour mémoire, la figure 1 rappelle de
l’évolution des paramètres en fonction des conditions locales et de
façon schématique l’allure des sections efficaces microscopiques
l’épuisement peut présenter des différences notables d’un code à
d’absorption de l’uranium 235 en fonciton de l’énergie.
l’autre. En revanche, l’intégration de ces codes dans l’environne-
ment plus général des châines de calcul peut se faire de façon très
variée en fonction des expériences et des exigences propres à
chaque concepteur (§ 7).
La figure 2 montre, de façon schématique, le découpage entre
les deux niveaux d’hétérogénéité dans un calcul type : celui d’un
cœur de REP comportant du combustible MOX, mélange d’oxydes
d’uranium et de plutonium. Elle met en évidence la structure très
hétérogène du crayon combustible, avec la succession de quatre
régions aux dimensions et aux propriétés différentes : la pastille, le
jeu interstitiel, la gaine et, extérieurement, le modérateur-réfrigé-
rant. La composition et la géométrie du crayon combustible déter-
minent l’importance de l’autoprotection des résonances (**). La
modélisation de ce phénomène, extrêmement important dans la
physique des réacteurs nucléaires, est une fonction majeure des
codes de cellule.
Figure 1 – Allure de la section efficace microscopique (**) Les discontinuités locales, brusques et parfois importantes, des propriétés neutro-
niques de la matière nucléaire sont dues, pour l’essentiel, à la nature résonnante des sec-
d’absorption de l’U235 (les résonances à très haute énergie ne sont tions efficaces microscopiques des isotopes et à la distribution spatiale des matériaux. La
pas représentées) prise en compte, fidèle et précise, de ces variations n’est possible actuellement que dans les
calculs probabilistes, qui admettent une description continue des propriétés énergétiques et
spatiales des milieux. Dans les calculs déterministes, ces inhomogénéités sont lissées, en
espace et en énergie, dans les processus de mise en groupes et d’homogénéisation des sec-
tions efficaces.
Le nombre moyen d’interactions des neutrons avec la matière est affecté par ces
approximations ; on constate alors une dérive des taux de réaction, variable avec le phé-

Toute reproduction sans autorisation du Centre français d’exploitation du droit de copie est strictement interdite.
© Techniques de l’Ingénieur, traité Génie nucléaire B 3 070 − 5
MÉTHODES DE CALCUL NEUTRONIQUE DE CŒUR _____________________________________________________________________________________________

Figure 2 – Deux niveaux d’hétérogénéité dans les calculs neutroniques : cas d’un REP à combustible MOX

nomène considéré. Pour contrebalancer cette tendance on adopte, dans les calculs, des certains domaines d’énergie. Les taux de fission et d’absorption
procédures spécifiques capables de rétablir, tout au moins partiellement, le bilan des
réactions neutroniques. Ces procédures se regroupent, dans l’essentiel, en deux grandes
neutroniques en sont modifiés d’autant (figure 3).
familles :
— celles qui adoptent localement un raffinement du maillage en énergie auquel on
superpose une description très détaillée des milieux physiques ;
— celles qui font appel au calcul de l’autoprotection des résonances. On utilise alors
des algorithmes appropriés qui ont pour objectif la recherche de sections efficaces
2.3 Génération des sections efficaces
effectives, obtenues grâce à l’application d’un principe d’équivalence. Ce formalisme, qui
s’avère particulièrement efficace pour le calcul des REP, comporte l’évaluation d’un ou plu-
pour les calculs de cœur
sieurs paramètres, dits facteurs d’autoprotection, dont le calcul peut être parfois
complexe.
L’autoprotection s’articule en deux volets, dont le principal est directement lié à la
composante énergétique, l’autre étant associé à la composante spatiale. Dans le cas d’un
Les sections efficaces microscopiques de tous les isotopes
barreau combustible d’UO2 , placé en milieu infini dans un environnement de type REP, présents dans le cœur et dans son environnement immédiat sont
l’autoprotection de l’uranium 238 s’élève à quelque 30 000 × 10–5, dont seulement 15 % compilées dans des bibliothèques, appelées aussi librairies, issues
sont dus à la composante spatiale. d’organismes internationaux. Actuellement, on fait référence à la
La figure 2 montre aussi, très schématiquement, la sensibilité de bibliothèque JEF-2 (Jointly Evaluated File ) livrée en 1992, résultat
la distribution des neutrons en énergie, dit spectre des neutrons, d’une collaboration entre les pays d’Europe et le Japon, élargie
aux caractéristiques des milieux ; ainsi la présence de forts absor- récemment à d’autres partenaires dont les États-Unis. Les données
bants, les grappes de contrôle par exemple, et la nature du contenues dans ces bibliothèques ne sont pas utilisables directe-
combustible, oxyde d’uranium ou oxyde mixte d’uranium et pluto- ment. Un traitement préalable, comportant l’intégration des para-
nium, modifient sensiblement la population neutronique dans mètres de résonance, la mise en groupes et la condensation en

Toute reproduction sans autorisation du Centre français d’exploitation du droit de copie est strictement interdite.
B 3 070 − 6 © Techniques de l’Ingénieur, traité Génie nucléaire
____________________________________________________________________________________________ MÉTHODES DE CALCUL NEUTRONIQUE DE CŒUR

efficaces en fonction d’un certain nombre de paramètres : entre

autres, la quantité d’énergie produite par le combustible (taux de
combustion), la présence de bore soluble et/ou de xénon, les tempé-
ratures du combustible et du modérateur et la densité du réfrigé-
rant.

2.4 Approximations numériques

La discrétisation d’un problème défini et borné sur un domaine

limité comporte toujours la définition d’un sous-ensemble de
points, segments ou sous-espaces appartenant à ce même domaine
auxquels on limite le calcul explicite des grandeurs recherchées
numériquement.
Les relations à établir entre les inconnues aux endroits choisis
(points de maillage) dépendent de la formulation mathématique
des problèmes à traiter et de leurs conditions aux limites. Elles ne
Figure 3 – Taux de fission en fonction de la léthargie sont jamais strictement exactes : le processus de discrétisation
introduit des erreurs de natures variées (troncature dans les séries,
arrondi dans les nombres, etc.) qui, tout particulièrement dans les
énergie, doit précéder la mise en forme des données requises par techniques aux différences finies, affectent la solution courante.
les codes de cellule (pour une information complète [3] [6]). Une fois la convergence atteinte, la solution obtenue s’écarte alors
En France, cette tâche est assurée par le couple de codes Njoy- d’autant plus de la solution exacte que les approximations liées à
Themis qui, dans le cas des REP, génèrent deux types de biblio- la discrétisation sont sévères. La solution obtenue par voie numé-
thèques, une à 99 groupes, une autre, plus détaillée, à 172 groupes rique ne rejoint la solution exacte que, à la limite, pour un maillage
d’énergie. Cette dernière bibliothèque possède la particularité d’être de calcul infiniment fin [5].
accessible par deux codes de cellule parmi les plus performants et Il est possible de pallier, au moins partiellement, les imprécisions
réputés en ce qui concerne le calcul des REP : Apollo-II pour la dues aux effets de maillage en couplant la discrétisation avec un
France, et Wims-E, pour le Royaume-Uni. développement polynominal des solutions recherchées entre les
Les bibliothèques utilisées par les codes de cellule sont calculées points (nodes ). De nombreux progrès ont été accomplis ces der-
une fois pour toutes et confinées dans des fichiers fournis en même nières années dans ce domaine grâce à l’évolution des techniques
temps que les codes de calcul. Il arrive toutefois que le concepteur de résolution. On dispose maintenant d’un large éventail de métho-
puisse réclamer une mise à jour ou une implémentation de ces des qui peuvent être regroupées en quatre familles principales : les
bibliothèques : la tâche est alors exécutée par les organismes différences finies, les éléments finis, les éléments finis mixtes et les
agréés (en France, le CEA). méthodes nodales.
Pour les RNR, le code Ecco, développé conjointement par les Par la suite (§ (4 et 5), nous reviendrons sur ces méthodes. Pour
Français et les Anglais dans le cadre des études EFR (European Fast l’instant, limitons-nous à quelques remarques sur les
Reactor ), fait appel à des bibliothèques très détaillées (plusieurs conséquences, en termes de précision et de représentativité, de la
milliers de micro-groupes). discrétisation de l’équation de la diffusion des neutrons.
L’analyse des spécificités des codes de cellule et de leur fonction- Nous avons vu ci-avant que toute discrétisation induit une impré-
nement sort du cadre de ce paragraphe. Nous nous limitons ici à cision variable avec le pas du maillage imposé au problème : en
quelques remarques et commentaires à caractère général, qui per- général, plus les mailles sont larges, plus la solution obtenue s’éloi-
mettent d’aborder une problématique vaste et complexe. Pour plus gne de celle que l’on obtiendrait avec un maillage infiniment fin. Par
d’informations sur ce sujet fondamental en physique des réacteurs, ailleurs, un rafinement très poussé du maillage n’est pas réalisable,
on renvoie à la nombreuse et très riche littérature internationale et, en pratique. Pour pallier cet inconvénient, les méthodes classiques
notamment, en ce qui concerne le code Apollo-II à [19] [20]. Pour aux différences finies utilisaient des techniques d’extrapolation
un exemple d’application d’Apollo-II aux calculs de projet dans un linéaire, par exemple celles développées par Kato [23]. Les métho-
contexte industriel, on renvoie à [21] et pour le code Ecco à [22]. des plus modernes, en revanche, font appel à des développements
polynomiaux ou en série de fonctions. Elles permettent ainsi
En général, les codes de cellule les plus modernes présentent un d’obtenir, avec des maillages relativement grossiers, des solutions
certain nombre de caractéristiques communes, à savoir : qui approchent de très près la solution asymptotique recherchée.
— ils résolvent l’équation du bilan neutronique dans sa forme
Ce n’est pas pour autant que l’on peut considérer le problème
intégrale en théorie du transport, généralement à deux dimen-
de la représentativité des résultats de calcul comme définitivement
sions. Des corrections sont prévues, s’il y a lieu, pour tenir compte
résolu. En effet, il reste toujours à vérifier comment les résultats se
de la dimension manquante ;
situeraient par rapport à ceux d’une expérience conceptuelle et/ou
— ils traitent des géométries hétérogènes complexes, avec un
à ceux d’un hypothétique calcul exact de référence qui reproduirait
nombre élevé de milieux aux propriétés physiques différentes, à
de façon fidèle la solution analytique, généralement impossible à
un nombre élevé de groupes d’énergie ;
calculer.
— ils comportent la possibilité de faire appel à plusieurs niveaux
d’approximation dans la résolution de l’équation du bilan, notam- À ce propos, il ne faut pas oublier que, dans le processus qui per-
ment en ce qui concerne le couplage des milieux et la probabilité met la génération des constantes utilisées dans les calculs de
d’interaction des neutrons avec la matière ; cœur (§ 2.3) à partir des données de base, sont effectués en
— ils possèdent des algorithmes sophistiqués pour le traitement cascade : plusieurs condensations en énergie, plusieurs homo-
des problèmes d’autoprotection des résonances, y compris la prise généisations en espace et un développement au premier ordre du
en compte de résonances multiples ; courant neutronique. Les approximations correspondantes ont un
— ils génèrent, si besoin est (c’est le cas pour les REP), des tables impact plus ou moins important sur les résultats selon qu’un pro-
à plusieurs niveaux de paramètrisation (tables multiparamétrées), cessus d’ajustement (équivalence) est mis en place ou non : la préci-
qui permettent de corriger par interpolation les valeurs des sections sion dépend étroitement des schémas adoptés et de la succession

Toute reproduction sans autorisation du Centre français d’exploitation du droit de copie est strictement interdite.
© Techniques de l’Ingénieur, traité Génie nucléaire B 3 070 − 7
MÉTHODES DE CALCUL NEUTRONIQUE DE CŒUR _____________________________________________________________________________________________

des étapes approximation-simplification-équivalence. Sous l’évolution du combustible et sa modélisation dans les codes de
certaines conditions, on peut, notamment, démontrer que, en raison calcul, se reporter à [3] [4] [5] ;
des compensations d’erreurs, le résultat d’un calcul de diffusion à — d’autre part, pour le calcul des REP, de fournir les informa-
mailles centrées, avec un seul point de calcul par région peut être, tions locales nécessaires à la définition des paramètres d’interpo-
même en présence de forts gradients, très proche de la solution lation. En raison du couplage étroit existant entre la production de
analytique de référence [24]. puissance et les conditions physiques des milieux dans les REP
(températures, densité du modérateur, etc.), ces informations sont
élaborées par des fonctions des codes, simulant un calcul thermo-
hydraulique sous forme simplifiée, ou bien elles sont transmises
2.5 Corrections ponctuelles directement par un module de calcul thermohydraulique couplé
des opérateurs avec le code de cœur.
Les données ainsi calculées permettent l’évaluation des para-
mètres physiques locaux utilisés, moyennant interpolation dans les
Nous avons vu que l’équation d’état d’un système nucléaire mul- tables multiparamètrées préparées par le code de cellule* pour
tiplicateur, après condensation, homogénéisation et discrétisation mettre à jour, à chaque pas de temps et en chaque node, les
des opérateurs de l’équation de Boltzmann sur un maillage donné, composantes des opérateurs définissant l’équation d’état du
peut s’écrire [équation (13)] : système multiplicateur.
A 0 Φ 0 + λ0 P 0 Φ 0 = 0 (*) La génération des tables multiparamétrées est effectuée en faisant varier de façon
couplée ou indépendante, toutes choses égales par ailleurs, un certain nombre de para-
mètres, choisis parmi les plus significatifs, autour du point de fonctionnement théorique
Cette équation comporte, outre les hypothèses réductrices sur la retenu pour le calcul nominal des constantes nucléaires, et en recalculant les sections effi-
forme de la solution que nous avons déjà discutées précédemment, caces dans ces conditions. Le nombre de paramètres considérés dépend du niveau
deux niveaux supplémentaires d’approximation qui touchent aux d’approximation choisi et des caractéristiques du code en aval. En général, pour les REP,
valeurs des composantes des opérateurs, à savoir : sont pris en compte, a minima, les températures du combustible et du modérateur, la
concentration de bore dans le modérateur et l’empoisonnement xénon. La génération de
— l’indépendance sous forme explicite du système par rapport ces tables est effectuée, de façon indépendante, à toute étape d’évolution.
au temps ;
— la connaissance parfaite et généralisée des composantes des
opérateurs dans tout le système. Ordres de grandeur des temps
et fraction de neutrons retardés en jeu
La première approximation implique que l’évolution du système
dans le temps, s’il y a lieu, est assimilée à une succession infinie
Temps de vie des neutrons prompts ......................... << 10–3 s
d’états stables, chacun très peu différent de son prédécesseur et
Valeur de référence pour les calculs
tout à fait indépendant et découplé de lui. En absence d’interven-
de conception des REP à combustible UO2 .............. 30 µs
tions externes qui viendraient perturber le système, cette approxi-
mation adiabatique est parfaitement justifiée en raison de l’écart Temps de vie des précurseurs de neutrons différés >1s
existant entre les constantes de temps des phénomènes en jeu Temps de vie des familles dont la vie
(pour plus d’information sur les cinétiques chimiques, se reporter est la plus longue......................................................... ≈ 1 min
au chapitre Théorie des réacteurs nucléaires, [B 3 025] et [2] [7]) :
— le retard avant l’apparition des neutrons de fission, quelques Retard avant l’apparition des effets
dizaines de secondes pour ceux dont l’attente est la plus longue de contre-réaction ........................................................ >1s
(moins d’un pourcent du total), une fraction de seconde (< 10–3 à Période des oscillations engendrées
10–7 s, selon le système considéré), pour tous les autres ; par les produits de fission à vie courte..................... > 10 h
— le changement des propriétés nucléaires des milieux,
consécutif à la production d’énergie (épuisement du combustible et Période de stabilité des propriétés neutroniques
accumulation des produits de fission) plusieurs heures, voire des matériaux par rapport
quelques jours avant que l’on puisse constater le moindre effet à l’épuisement du combustible nucléaire ................. ≈ 1 mois
significatif.
Fraction de neutrons retardés :
La seconde approximation concerne, pour chaque palier, les opé-
dans les REP (en moyenne) ........... entre 450 et 700 × 10–5
rateurs qui décrivent le système et qui sont supposés connus et
dans les RNR.................................... de l’ordre de 350 × 10–5
déterminés en tout point.
Si cette connaissance est parfaitement envisageable d’un point
de vue théorique, dans la pratique, en raison de considérations de
coût, d’encombrement mémoire et d’efficacité de calcul, seules les
sections efficaces d’un nombre très restreint de points sont acces- 3. Calculs du transport
sibles directement. Toutes les autres peuvent être reconstruites
grâce à des algorithmes qui prennent en compte les valeurs locales des neutrons
des paramètres. Les méthodes plus modernes et plus sophisti-
quées font appel à l’interpolation dans les tables multiparamétrées,
d’autres à des développements analytiques, d’autres encore à une
combinaison des deux. 3.1 Méthodes déterministes (analytiques)
Les deux approximations décrites ci-dessus nécessitent la pré-
sence dans les codes de calcul de cœur de sous-programmes Dans la plupart des cas, il n’est pas possible de résoudre direc-
capables : tement par voie analytique l’équation générale du transport (cf.
— d’une part, de mettre à jour les données relatives aux Annexe), mais on peut en rechercher des solutions approchées par
compositions isotopiques du combustible nucléaire afin d’y inté- diverses méthodes numériques.
grer les effets de la production d’énergie. Il s’agit là, en général, de
En faisant abstraction de la variable temps, le flux est une fonc-
sous-programmes dits de calcul d’évolution qui, à l’aide de formu-
tion de six variables : une pour l’énergie, deux pour la direction
lations analytiques, ou, plus souvent, de méthodes de résolution
angulaire et trois pour l’espace. Le plus souvent, on cherche une
matricielles, résolvent, par node, les équations d’évolution des
solution indépendante du temps (on ne traite explicitement ce der-
nuclides sous une forme simplifiée ; pour plus d’information sur

Toute reproduction sans autorisation du Centre français d’exploitation du droit de copie est strictement interdite.
B 3 070 − 8 © Techniques de l’Ingénieur, traité Génie nucléaire
____________________________________________________________________________________________ MÉTHODES DE CALCUL NEUTRONIQUE DE CŒUR

nier que pour l’analyse de situations transitoires), ce qui revient à une somme finie de termes du développement (formule de quadra-
déterminer le flux Φ (r, E, Ω) en tout point d’un milieu de volume V. ture).
■ Discrétisation par rapport à l’énergie, calcul multigroupe Ces polynômes sont orthogonaux : 〈 P œ ⋅ P m〉 = 0 si œ ≠ m , ce qui
On discrétise la plage d’énergie dans laquelle se trouvent les conduit à exprimer les sections efficaces et le flux suivant le
neutrons (de plusieurs millions d’électronvolts, partie haute du formalisme :
spectre de fission, à quelques millièmes d’électronvolt, partie L

∑  ----------------
2œ + 1 œ
basse du spectre) en intervalles E g B E B E g + 1 (§ 2.2). On réduit Σ(r, E′ → E, µ) = Σ ( r , E ′ → E ) Pœ ( µ ) (19)
4π 
ainsi un problème dépendant de l’énergie en un système d’équa- œ=0


+1
tions couplées : équations multigroupes, chacune d’entre elles œ
ayant la forme d’une équation mono-énergétique, où apparaissent avec Σ ( r, E′ → E ) = 2π Σ ( r, E′ → E, µ ) P œ ( µ ) d µ
des valeurs du flux Φg (r, Ω) et des sections efficaces Σg (r, Ω) par –1
L
groupes d’énergie :
∑  ----------------
2œ + 1
Φ ( r, E, µ ) = Φ œ ( r, E ) P œ ( µ )

et (20)
Eg + 1
4π 
Φ g ( r, Ω ) = Φ ( r, E, Ω )dE œ=0


Eg +1
(16)
avec Φ œ ( r, E ) = 2π Φ ( r, E, µ ) P œ ( µ ) d µ

Eg + 1
Σ ( r ,E ) Φ ( r, E, Ω ) –1
et Σ g ( r, Ω ) = ----------------------------------------------dE
Φ g ( r, E, Ω )
Eg En règle générale, les sources sont traitées comme étant isotro-
pes.
■ Discrétisation par rapport à l’espace
En portant les expressions (19) et (20) dans l’équation générale
On découpe le volume à étudier en un ensemble de volumes élé- du transport, on obtient un système de N + 1 équations différentiel-
mentaires (maillage) qui peuvent chacun être considérés comme les.
quasi homogènes et, suivant les cas, on calcule un flux au centre
ou à la périphérie des mailles. Comme nous l’avons vu précédem- Dans la solution globale, l’importance du polynôme diminue rapi-
ment, le choix du maillage est lié à la technique de résolution dement lorsque son degré augmente, aussi limite-t-on le dévelop-
numérique que l’on met en œuvre. pement aux premiers termes. Si l’on s’arrête à  = N , on dit que
l’on utilise l’approximation PN ; la pratique conduit à choisir N
■ Discrétisation par rapport à l’angle impair.
La discrétisation du problème par rapport aux directions angu- On peut montrer que limiter le développement aux harmoniques
laires peut se faire : d’ordre 0 et 1 (approximation P1) ramène à la théorie de la diffusion
— en tenant compte de la distribution des vitesses des neutrons (cf. chapitre Ralentissement et diffusion des neutrons [B 3 014]).
après leur naissance et/ou à l’issue des chocs ; dans ce cas, on
■ Méthodes des ordonnées discrètes
développe les sections efficaces de ralentissement avec le flux en
polynômes et on découpe l’espace des directions en secteurs : Dans cette méthode [25], on remplace les directions Ω’ et Ω par
c’est l’approche intégro-différentielle ; un ensemble limité de directions Ωi .
— en considérant seulement la probabilité que les neutrons, en se Soit Φi (r, E, Ωi ), le flux angulaire dans la direction Ωi , on va
propageant en ligne droite, puissent atteindre un point quelconque choisir les directions Ωi , leur nombre N et leur poids respectifs ωi
du système avant d’interagir avec la matière : c’est l’approche inté- pour que :


grale. N
Φ ( r, E, Ω ) dΩ = ∑ ωi Φ i ( r, E, Ω i ) (21)
4π i=1
3.1.1 Méthodes intégro-différentielles
Quand le flux présente une symétrie azimutale, on repère les Ωi
■ Développement sur des harmoniques sphériques par le seul paramètre Θ et l’on utilise µi = cos Θi , avec µ0 = – 1 et
µN = + 1. On choisit pour valeur de µi les N racines du polynôme de
Dans l’espace, on repère la direction Ω par deux paramètres, les
Legendre de degré n (n pair) et pour ωi les solutions de :
angles θ et ϕ : lorsque le système présente une symétrie azimutale,
comme dans le cas des géométries plane ou sphérique, un seul N
2
∑ ωi
angle suffit à repérer la direction. n
µ i = ------------- avec n = 0, 2, 4 ... (2N – 2) (22)
n+1
Une méthode couramment utilisée consiste à développer le flux i=1
angulaire et les sections efficaces sur une base de fonctions appro-
priées, en général des harmoniques sphériques suivants : Dans les méthodes intégro-différentielles, on fait l’approximation
d’une variation linéaire de Φ (r, E, Ω) entre les Ωi ; on peut ainsi expri-
( 2œ + 1 ) ( 1 – m ) ! m im ϕ
mer simplement les termes de l’équation intégro-différentielle fai-
Y œ ,m ( θ , ϕ ) = ----------------------------------------------- P œ ( cos θ ) e (17) sant intervenir les dérivées par rapport aux coordonnées
4π ( 1 + m ) !
angulaires [26]. La précision et la complexité des calculs aug-
avec œ et m entiers. mentent avec N ; en général, on se limite à N = 4 ou 6, méthodes
appelées S4 ou S6 , parfois à N = 8, méthode S8 , rarement au-delà.
Si le milieu peut être considéré comme isotrope, ce qui est le plus
Si on fait l’hypothèse que le flux est linéaire à la fois en fonction
souvent le cas dans les calculs de cœur, alors la section efficace
de l’espace et de l’angle sur chaque intervalle de la discrétisation,
Σ (r, E’ → E, Ω’ → Ω) ne dépend que d’un paramètre : l’angle Θ que
on peut mettre en œuvre une technique de résolution numérique
font Ω’ (après choc) et Ω (avant choc), et que l’on va repérer par
basée sur le schéma diamant [3] [4].
µ = cos Θ :
2 3
3 µ –1 5µ – 3µ
P0 (µ ) = 1, P1 (µ ) = µ, P2 (µ ) = ----------------- , P 3 ( µ ) = ----------------------- (18)
2 2 3.1.2 Méthodes intégrales (probabilité de collision)
Les P œm sont les fonctions associées de Legendre qui pour m = 0 Si l’on suppose les sections efficaces indépendantes du temps,
se réduisent aux polynômes de Legendre P œ ( µ ) . On peut alors alors les neutrons qui arrivent en r, à l’énergie E et dans la direction
remplacer dans l’équation intégro-différentielle les intégrales par

Toute reproduction sans autorisation du Centre français d’exploitation du droit de copie est strictement interdite.
© Techniques de l’Ingénieur, traité Génie nucléaire B 3 070 − 9
MÉTHODES DE CALCUL NEUTRONIQUE DE CŒUR _____________________________________________________________________________________________

Ω ont été créés dans l’ensemble des points r’ tels que r’ = r – x Ω, Toutefois, leur utilisation demande d’éviter l’apparition de phé-
et ils ont parcouru la distance r’ – r sans subir de choc. nomènes numériques tels que la distorsion du flux calculé dans
On passe de la forme intégro-différentielle à une forme intégrale certaines directions privilégiées par le calcul (effet de raies).
en introduisant le terme : Les méthodes des probabilités de collision sont utilisées pour la

Σ
r détermination des constantes neutroniques utilisées dans les cal-
Σ̃ tx = – t ( r′, E, Ω )dx culs de cœur (§ 2.2), ce qui se fait avec des codes de réseaux tels
A que Apollo [32] [33] en France ou Wims [34] en Grande-Bretagne.
probabilité de propagation en ligne droite entre le point de départ
A et le point d’arrivée générique r.
On montre alors [3] que l’équation du transport se met, dans le
3.2 Méthodes probabilistes (Monte-Carlo)
cas d’un milieu à une dimension, sous la forme :


∞ Dans le paragraphe 3.1, on a vu que les erreurs systématiques
˜
Φ ( r, E, Ω ) = e– Σ tx q ( r′, E, Ω ) dx (23) provenaient non seulement des incertitudes sur les données fonda-
0 mentales de base (sections microscopiques) mais aussi des
approximations méthodologiques. Parmi ces dernières, on peut
avec q(r’, E, Ω) représentant le nombre de neutrons d’énergie E et citer la séparation espace-direction angulaire-énergie des neutrons
de direction Ω émis en r’ = r – x Ω, auquel on ajoutera éventuelle- dans les calculs spectraux, la non-utilisation systématique de modé-
ment les sources. lisation tridimensionnelle et l’approximation énergétique multi-
Ce formalisme se généralise à plusieurs dimensions. groupe dans les études de projets industriels. Ces études sont
Dans la méthode de résolution par probabilité de collision [27], relatives à des calculs de cœurs de réacteurs, des piscines de
on découpe l’espace en éléments de volumes V i , que l’on va stockage d’assemblages usés, des calculs de protection et toutes
considérer comme homogènes, et l’on calcule les facteurs de pro- autres études mettant en jeu des neutrons.
babilité de collision ij qui représentent la probabilité qu’un neu- À l’opposé, les méthodes de Monte-Carlo [8] sont capables de
tron émis dans le volume Vi subisse sa première interaction dans simuler des configurations complexes en trois dimensions en trai-
le volumeVj . On entend par neutron émis dans Vi tout neutron qui tant de façon continue et simultanée l’énergie, l’angle et la direction
y a été créé ou dont les caractéristiques (énergie, direction) y ont des neutrons. Ces méthodes permettent de lever certaines incerti-
été modifiées suite à une diffusion ou à une autre interaction. tudes propres aux méthodes déterministes. Les erreurs résultant
On décrit ci-dessous le principe de la méthode dans le cas où des méthodes sont bien sûr celles implicites aux données fonda-
l’on suppose les chocs isotropes. mentales de base ainsi que celles inhérentes à la méthode qui pren-
nent la forme d’une erreur statistique. En effet, comme son nom
On note Σ t et Σ s les sections efficaces totale et de diffusion
l’indique, la méthode de Monte-Carlo est une méthode de résolution
dans le volume générique V,
numérique basée sur l’application de propriétés des variables aléa-
P et Φ la source moyenne et le flux moyen dans V. toires et sur l’utilisation de techniques d’échantillonnage au hasard.
On commence par calculer les probabilités de collision ij : Elle fait appel aux méthodes statistiques et aux probabilités.

vj
 3
d r′ Σ t ( r′ ) 
vj
3 e tx
– Σ̃
d r ------------2- Σ s ( r ) Φ ( r ) + P ( r )
4πx
3.2.1 Principe
 ij = --------------------------------------------------------------------------------------------------------------------------- (24)
 vj
3
d r Σs ( r ) Φ ( r ) + P ( r )
La mise en œuvre de la méthode commence par la construction
d’histoires de particules en considérant la succession des événe-
ments suivants :
Si Vi est petit, on pourra supposer que la densité des départs de — émission de particules sources (répartition spatiale, énergétique
neutrons y est uniforme d’où : et angulaire) ;
— choix d’un parcours dans un milieu de section efficace
Σ tj 3
 ij = ------ d r ′
Vi vi
 vi
3 e tx
– Σ̃
d r -------------2-
4π x
(25)
totale Σ t ;
— nature du noyau rencontré et type d’interaction ;
— diffusion isotrope ou anisotrope.
Puis, on calcule le flux moyen dans chacun des volumes par : L’histoire de la particule est poursuivie jusqu’à sa disparition soit
par absorption dans le système, soit par fuite. Chaque événement
Vi Σ ti Φ i = ∑ ij Vi ( Σ sj Φ j + Pj ) (26) décrit ci-dessus est traité d’une manière aléatoire en utilisant des
nombres au hasard (tirage) et en comptabilisant certaines quantités
j
numériques relatives aux situations rencontrées au cours des his-
Ce calcul se résume à la résolution d’un système linéaire de la toires ainsi construites. C’est donc le dépouillement statistique de
forme (26). ces histoires qui donne les informations recherchées [9].
Il est possible de prendre en compte une anisotropie linéaire des ■ Tirage selon une loi de probabilité quelconque
chocs en modifiant la formulation (24) par une correction de
En général, les événements aléatoires étudiés ne correspondent
transport [28].
pas à une distribution équiprobable mais plutôt à une distribution
quelconque ƒ (x ) sur un domaine [a,b] de x. Par définition :
3.1.3 Domaines d’application

b
ƒ (x ) d x = 1 (27)
Les méthodes intégro-différentielles sont particulièrement utili- a
sées pour le calcul du flux dans des milieux hétérogènes où les
variations spectrales du flux sont importantes. Ces méthodes sont On peut tout aussi bien considérer l’intégrale :


relativement simples à mettre en œuvre dans le cas des géomé- x
tries à une dimension, code ANISN [29] ; elles deviennent plus ƒ(x’ ) dx’ = F (x) (28)
complexes pour les cas à deux ou trois dimensions, codes Dot et 0
Tort [30] [31].

Toute reproduction sans autorisation du Centre français d’exploitation du droit de copie est strictement interdite.
B 3 070 − 10 © Techniques de l’Ingénieur, traité Génie nucléaire
____________________________________________________________________________________________ MÉTHODES DE CALCUL NEUTRONIQUE DE CŒUR

qui exprime la probabilité que la variable aléatoire prenne une <g > et <g 2> ne sont pas connus et donc σ N ne peut être calculé.
valeur inférieure ou égale à x, F (x) varie toujours de 0 à 1. On peut néanmoins l’estimer en posant :
La méthode de Monte-Carlo repose sur le fait qu’à partir d’une
2 2
séquence de nombres équiprobables sur [0, 1] il est possible <g > = g
d’atteindre une variable aléatoire distribuée selon ƒ(x ) sur [a, b]. Il
suffit de tirer des nombres si selon une séquence équiprobable et 2 N 2 1 2
<g> = -------------- ( g ) – ------ g (34)
de calculer xi correspondant par l’équation : N–1 N
si = F (xi ) (29) 2
(à noter que ( g ) est un estimateur biaisé de <g 2>), il vient :
Exemple : soit : 2
2
g – ( g)
ƒ(x ) dx = K exp (– Kx ) dx ( σN )2 ≈ ------------------------
N–1
- (35)
avec K constante de proportionnalité
On voit ainsi que l’erreur probable commise en prenant l’estima-
La probabilité pour un neutron parti de x = 0 de subir une réaction tion g N à la place de <g > varie comme 1 ⁄ N (racine carrée de
entre x et x + dx. la variance ou écart-type). Pour réduire cette erreur d’un facteur 10,
On a : il faut étudier 100 fois plus d’histoires. C’est un inconvénient


x sérieux de la méthode de Monte-Carlo, aussi l’a-t-on améliorée en
F(x) = ƒ(x’ ) dx’ = 1 – exp (– Kx) utilisant des techniques de réduction de variance.
0

Si on tire si , on calcule :
1
3.2.2 Techniques de réduction de la variance
x i = – ----
-
K ln ( 1 – s i )
Pour que l’écart-type soit faible, il faut donc effectuer un très
x i est la distance à laquelle a lieu l’événement réaction dans l’histoire grand nombre de tirages, ce qui se traduit par un grand nombre de
de ce neutron. particules devant atteindre et, éventuellement, traverser la cible
d’intérêt. Afin de réduire le nombre de tirages (histoires), on utilise
On procède ensuite, par exemple, au choix entre une réaction de des techniques particulières dites techniques de réduction de la
diffusion ou une réaction d’absorption en tirant un nouveau s. Les variance. Deux des techniques les plus souvent utilisées sont le
sections étant Σa et Σs avec Σ = Σa + Σs , si s B Σa /Σ on retient qu’il splitting [35] [36] et la pondération par importance.
y a absorption sinon qu’il y a diffusion. On poursuit l’analyse en
recherchant la direction du nouveau parcours, puis la distance par- Il faut noter que, si la méthode de Monte-Carlo est utilisée pour
courue avant la nouvelle réaction et ainsi de suite. On enregistre au obtenir une réponse en réactivité (milieu multiplicateur de neutrons)
fur et à mesure toutes les quantités pour lesquelles on désire des d’un système quelconque, tous les neutrons doivent participer au
résultats. processus. Pour minimiser l’écart-type, la seule solution est donc
le traitement d’un nombre suffisant et, par conséquent, important
■ Analyse statistique d’histoires (quelques millions). Les techniques de réduction de
Si l’on veut calculer la vraie moyenne : variance qui sont discutées ci-après ne sont donc appliquées qu’aux
problèmes spécifiques de protection ou de propagation de parti-

 cules hors du cœur du réacteur.

+∞
<g> = g (x ) ƒ(x ) dx (30)
–∞
■ Splitting
Dans les calculs de propagation profonde, les effets d’angle
On tirera x selon ƒ(x ) pour obtenir des xi , avec lesquels on cal- solide font en sorte que la plupart des particules ratent la cible
cule la moyenne : (qui, en général, possède des dimensions assez petites par rapport
N à celles de la source et aux distances en jeu) : les temps de calcul
1
g
N
= ------
N ∑ g ( xi ) après N tirages (31) risquent alors de devenir prohibitifs. La technique du splitting, qui
consiste à dédoubler en fragments de poids réduit chaque parti-
i=1
cule franchissant une surface donnée dans la direction de la cible,
Pour N fini, g N (quantité aléatoire) est en général différent de permet de conserver partiellement l’information de la composante
<g >. Le théorème de la limite centrale indique que la probabilité, de vitesse utile. Après génération, chaque fragment est suivi indi-
notée  { g N – <g> < ε }, que g N – <g> soit inférieur à ε tend, viduellement et peut, à son tour, subir des fragmentations.
pour N grand, vers l’expression : L’idée est donc de multiplier le nombre de particules, quitte à en

 ε N ⁄σ 2 diminuer le poids, aux environs de la cible. Chaque fragment tra-

1⁄2
2--- exp  – ----- d t
t
(32) versant la cible fournira, alors, une contribution — réduite — à la
 π 0
 2 fonction réponse recherchée. Le procédé complémentaire, la roulette
russe, est appliqué aux échantillons, ou fragments d’échantillons,
avec la variance σ 2 = <g 2> – <g >2. traversant une surface de splitting dans la mauvaise direction par
rapport à la cible : une fois sur deux, on annihile la particule et l’autre
Cela signifie que, reprenant plusieurs fois le calcul qui donne g N , fois on continue à la suivre avec un poids multiplié par deux. On
on obtiendrait autant de valeurs g N différentes mais distribuées génère ainsi une population éparse de particules lourdes loin de la
cible.
autour de la vraie valeur <g > selon une distribution normale. On
2 Les surfaces de splitting découpent la géométrie de la
montre que la variance ( σ N ) de gN vaut : configuration à étudier en cellules à chacune desquelles est attri-
buée une importance. Cette importance est la contribution attendue
2 2 2
σ d’une particule apparaissant dans cette cellule à la réponse d’un
( σN )2 ≈ -----
N
<g > – <g>
- ≈ ----------------------------------
N
(33) détecteur dans la cible. Mathématiquement, c’est la solution du pro-
blème adjoint dont la source est la fonction réponse du détecteur
dans la cible.

Toute reproduction sans autorisation du Centre français d’exploitation du droit de copie est strictement interdite.
© Techniques de l’Ingénieur, traité Génie nucléaire B 3 070 − 11
MÉTHODES DE CALCUL NEUTRONIQUE DE CŒUR _____________________________________________________________________________________________

■ Pondération par importance De nombreuses méthodes ont été développées pour atteindre
L’autre procédé de réduction de variance n’est pas fondamenta- cet objectif. Elles évoluent très rapidement et on assiste à l’heure
lement différent du splitting, il repose aussi sur l’utilisation de actuelle, avec le développement des ordinateurs, à une montée en
l’importance. puissance des modules résolvant directement l’équation du trans-
port dans des géométries complexes. Ces méthodes ne sont pour
La notion de pondération, c’est-à-dire la théorie des jeux biaisés l’instant pas parvenues à maturité et les concepteurs continuent à
en mathématique statistique, correspond dans notre cas à un flux faire appel principalement aux outils résolvant l’équation de la dif-
de particules biaisé : fusion neutronique.
Φ I ( x ) = Φ (x ) I ( x ) (36)
Comme il a été dit au paragraphe 2, les calculs de distribution
où la fonction I (x ), appelée importance, est choisie grande dans spatiale de neutrons à l’intérieur du réacteur comportent la traduc-
les régions de l’espace où sont recherchés les résultats. D’où tion de l’équation de la diffusion en un système d’équations algé-
l’expression suivante des flux biaisés : briques. En pratique, le domaine sur lequel l’équation d’état du
réacteur est définie et bornée est subdivisé spatialement en un
1 ensemble de régions dites mailles (ou points de calcul) dans les-
Φ I ( x ) = Φ ( x ) ----- exp ( Kx ) (37)
A quelles l’équation est intégrée spatialement. Ce procédé, appelé
discrétisation, permet de passer d’un problème aux dérivées par-
La notion d’importance est fondamentale et l’intérêt du biaisage tielles (dont la solution analytique n’existe pas pour un réacteur
est de réduire la variance sans altérer l’espérance mathématique, réel) à un problème numérique.
c’est-à-dire, dans le cas de la résolution du transport par une
méthode de Monte-Carlo, de traiter uniquement les particules Plusieurs techniques de discrétisation sont possibles.
contribuant au résultat recherché en éliminant ainsi la plus grande ■ L’une d’entre elles, utilisée encore actuellement pour les études
partie des particules inutiles. neutroniques de conception des réacteurs de puissance, est celle
Cette méthode demande des réglages des différentes impor- des différences finies (cf. chapitre Approximation des équations aux
tances données suivant les chemins parcourus par les particules. Le dérivées partielles. Méthodes aux différences finies [A 550] dans le
choix de l’importance n’est pas prédominant sur le résultat. Un traité Sciences fondamentales). Elle consiste à évaluer le courant de
mauvais choix donne un résultat correct mais au détriment d’un neutrons (exprimé par la loi de Fick) entre une maille et l’autre, dans
écart type plus important et/ou d’un temps de calcul détérioré. une direction donnée, de la façon suivante :
(Φ – Φ )
3.2.3 Domaines d’application J = – D grad ( Φ ) ≈ – D ------------------------
a b
- (38)
h
L’utilisation de méthodes de Monte-Carlo couvre le domaine où D est la constante de diffusion, a et b sont deux points du
complet allant du traitement des sections efficaces de base (neu- domaine où le flux Φ doit être évalué et h leur distance. Dans cette
trons, gamma, etc.) jusque, et cela grâce à la puissance croissante approximation, la dérivée spatiale est approximée, au premier
des ordinateurs, aux calculs de cœur ou de piscines de désactiva- ordre, par un rapport de différences. On présente plus de détails sur
tion. cette approximation dans le paragraphe 4.2.
Ces méthodes permettent de prendre en compte exactement les ■ Une autre technique utilisée en neutronique est celle des élé-
particularités d’une géométrie donnée et des interactions particules/ ments finis. Elle consiste à exprimer le flux comme développement
matière. Ces dernières sont : de fonctions de base (normalement des polynômes) :
— pour les neutrons : choc élastique avec anisotropie quelconque,
réactions inélastiques et réactions (n → 2n ), capture, fission et ther- Φ = a0 + a1 f1 (x, y, z ) + a2 f2 (x, y, z ) + ... (39)
malisation dans un gaz libre ou avec liaisons moléculaires ; Dans ce cas, la dérivée spatiale qui exprime le courant est for-
— pour les gamma : diffusion Compton, effet photoélectrique et mulée analytiquement. L’intégration dans chaque maille (avec pro-
production de paires. jection sur les fonctions orthogonales) conduit à des équations
L’utilisation des méthodologies Monte-Carlo comme référence dont les inconnues sont les coefficients de ces fonctions. Étant
pour la validation peut s’appliquer, entre autres, à des études de donné que cette technique n’est pas utilisée industriellement pour
problèmes neutrons ou gamma à courtes distances des sources les études de conception, elle ne sera pas détaillée ici. Pour plus
(échauffement, dommage, irradiation, etc.), de problèmes de péné- d’information, se reporter au chapitre Méthode des éléments finis
tration profonde en neutrons ou gamma dans des géométries [A 656] dans le traité Sciences fondamentales et [5] [10] [40] [41]
complexes et de problèmes critiques avec multiplication des neu- [42].
trons de fission. Citons quelques noms de codes les plus généra- La technique de discrétisation utilisée de préférence actuelle-
lement utilisés dans le monde : Tripoli en France [37] MCNP [38] et ment pour les études des REP est l’expansion nodale (§ 5). Pour les
Vim aux États-Unis [39]. RNR, on adopte toujours la technique des différences finies (§ 7.3).

4. Méthodes de résolution 4.2 Discrétisation et effet de maillage

de l’équation de la diffusion
La particularité de la technique de discrétisation aux différences
finies est d’approximer la dérivée spatiale du flux par un rapport
de différences du flux et de la distance. Cette approximation peut
4.1 Méthodes aux différences être effectuée de plusieurs façons [11] [12].
et aux éléments finis La plus immédiate consiste à exprimer le courant entre une
maille a et sa voisine b par l’intermédiaire d’un flux de surface Φs .
Les codes de cœur ont pour objectif principal l’évaluation de la Si on fait l’hypothèse que la valeur moyenne du flux dans une
réactivité des systèmes multiplicateurs et le calcul des distribu-
tions fines de puissance (et éventuellement, d’autres taux de réac-
tion crayon par crayon).

Toute reproduction sans autorisation du Centre français d’exploitation du droit de copie est strictement interdite.
B 3 070 − 12 © Techniques de l’Ingénieur, traité Génie nucléaire
____________________________________________________________________________________________ MÉTHODES DE CALCUL NEUTRONIQUE DE CŒUR

maille est égale à la valeur du flux en son centre, le courant à la Le choix de la taille de chaque maille conditionne la précision du
surface d’une région peut être écrit : résultat. Un maillage trop large ne permet pas de suivre des varia-
tions rapides de flux. Un maillage trop serré augmente par ailleurs
Da ( Φs – Φa ) Db ( Φb – Φs ) l’influence des erreurs de troncature sur le résultat, étant donné que
J = – --------------------------------- = – ---------------------------------- (40)
ha ⁄ 2 hb ⁄ 2 la différence de flux qui apparaît dans la relation (47) peut porter sur
les derniers chiffres significatifs*. En général, l’augmentation du
En éliminant le flux surfacique, on arrive à exprimer le courant nombre de mailles introduit des instabilités numériques qui aug-
par l’intermédiaire des flux dans les deux mailles : mentent le nombre d’itérations.
(*) L’effet de maillage est analysé en détail dans le paragraphe 2. Il est important
J = Rab (Φb – Φa ) d’observer qu’il n’est pas le seul effet de modèle qui doit être repris : l’hétérogénéité
énergétique et spatiale sont autant de sources d’incohérence des résultats des calculs de
avec : cœur.
Da Db
--------
- --------
ha hb
R ab = 2 ---------------------------- (41)
Da Db
- + --------
--------
4.3 Processus de calcul
ha hb

Une autre façon d’approximer la dérivée du flux pour représen- Le système homogène décrit dans le paragraphe 4.2 peut être
ter la loi de Fick consiste à écrire : écrit en forme compacte :
1
D m ( Φa – Φb ) A 0 Φ 0 = --------- P 0 Φ 0 (48)
J = – ------------------------------------- (42) k eff
ha hb
------ + ------ où A0 est l’opérateur matriciel qui décrit la diffusion, l’absorption
2 2
et le ralentissement des neutrons (il contient les termes Rij , Σa , Σr)
où Dm est une constante de diffusion pondérée selon la relation : et P0 est l’opérateur matriciel de production de neutrons (il contient
les termes de source).
h a Da + h b Db
D m = ---------------------------------
- (43) Ce système ne peut être résolu directement en raison de la taille
ha + hb des matrices A0 et P0.
la constante de couplage entre les mailles a et b étant alors : La méthode couramment utilisée consiste à passer par l’inter-
médiaire d’un système non homogène plus facile à résoudre. En
ha Da + hb Db effet, la relation (48) peut être séparée en deux équations, à savoir :
R ab = 2 ---------------------------------
2
- (44)
( ha + hb ) A 0 Φ 0 = P˜ (49)
Pour une maille à la frontière du domaine, le courant est formulé
1
à partir des conditions aux limites (§ 2.1). On considère ici la seule P˜ = --------- P 0 Φ 0 (50)
condition d’annulation du flux. Que ce soit l’une ou l’autre des k eff
façons d’approximer la dérivée du flux, le courant s’exprime tou-
Le processus itératif consiste alors à résoudre plusieurs fois le
jours par (Φs = 0) :
système (48) en réévaluant, à chaque fois, le terme de source P˜ à
Da Φa l’aide de l’équation (50).
J = – -------------- (45)
ha ⁄ 2 La valeur propre est calculée à chaque itération. Plusieurs métho-
des de calcul existent, parmi lesquelles on retiendra les suivantes :
et la constante de couplage devient : (n)
(n) 〈 P0 , Φ 〉
Da k eff = ---------------------------
(n)
(51)
R ab = ------
- (46) 〈 A0 , Φ 〉
ha
ou :
En introduisant la relation (42) dans l’équation homogène de la
(n)
diffusion (§ 2.1) et en intégrant sur le volume de la maille, le cou- (n) 〈 f, P 0 Φ 〉
rant disparaît. La résolution du problème est alors ramenée à celle k eff = -------------------------------
(n)
- (52)
d’un système de NG équations (N étant le nombre de mailles et G 〈 f, A 0 Φ 〉
le nombre de groupes de neutrons) où les inconnues sont les
valeurs du flux dans chaque maille pour chaque groupe de neu- ou encore :
(n)
trons. (n) 〈 P0 , Φ 〉 (n – 1)
k eff = ---------------------------------
(n – 1)
- k eff (53)
Dans ces conditions, en utilisant les définitions (8), l’équation (14) 〈 P0 , Φ 〉
devient, pour chaque maille i et groupe g :
L Le symbole <,> indique classiquement une intégration en espace
∑  Rgiœ  Φ g – Φ g  S iœ +  Σ ga, i + ∑ Σ gi → g′ Φ g0i V i et en énergie sur tout le domaine et n représente l’indice d’itéra-
 0i 0œ   tion.
œ=1 g′
La relation (51) est obtenue en intégrant l’équation (48) sur
∑  Σi + --------- χ ν Σ f ,i Φ 0i V i = 0
g′ → g 1 g g g′
– (47) l’espace et l’énergie. La relation (52) est obtenue en effectuant une
k eff 
g′ convolution de l’équation (48) sur une fonction générique de
g l’espace et de l’énergie. Cette fonction peut être redéfinie à chaque
où R iœ est la constante de couplage entre les mailles i et œ , définie
itération. Dans certains codes de calcul à deux groupes d’énergie,
par une des relations (41), (44) ou (46) et S iœ est la surface de
elle est supposée égale aux flux rapides (hypothèse dite du flux
contact entre ces deux mailles (L vaut 6 en géométrie rectangulaire
rapide autoadjoint ) et ne dépend pas de l’énergie. La relation (53)
3D).
est obtenue en substituant l’équation (50) relative à l’itération n – 1
Ce système, en général, ne peut pas être résolu directement : sa dans l’équation (48) relative à l’itération n.
solution est donc recherchée par une approche itérative.

Toute reproduction sans autorisation du Centre français d’exploitation du droit de copie est strictement interdite.
© Techniques de l’Ingénieur, traité Génie nucléaire B 3 070 − 13
MÉTHODES DE CALCUL NEUTRONIQUE DE CŒUR _____________________________________________________________________________________________

Avant la première itération, le flux est initialisé à une distribution d’appeler itération de diffusion l’itération interne (parce qu’elle
arbitraire par maille et par groupe. représente la diffusion des neutrons produits par la source P˜ ) et
Le processus itératif est arrêté lorsque l’erreur maximale en itération de source l’itération externe.
espace et en énergie est inférieure à un seuil préétabli ou défini par Il n’est pas nécessaire de faire converger totalement le processus
l’utilisateur (critère de convergence) : de diffusion à chaque itération de source. La procédure courante
est d’effectuer un nombre fixé par avance d’itérations de diffusion
(n) (n – 1)
Φi – Φi pour chaque itération de source et de ne tester que la convergence
(n – 1)
- <ε
max i ------------------------------------- (54) du processus externe.
Φi
La vitesse de convergence peut être augmentée par des métho-
Cet algorithme qui converge naturellement (dans les situations des d’accélération [11] [12] [13]. Notamment, pour réduire le nom-
courantes en neutronique) est connu sous le nom de méthode des bre d’itérations de diffusion, on applique la méthode de
puissances (cf. chapitre Méthodes numériques de base [A 1 220] surrelaxation qui consiste à appliquer une extrapolation aux flux
dans le traité Sciences fondamentales et [43]. calculés :
La résolution du système (49) et (50), bien que simple, n’est pas Φext = ω Φ (n ) + (1 – ω ) Φ (n – 1) (58)
possible en général de façon directe à cause de la taille importante
des matrices. Sa solution est obtenue par un autre processus itératif. où ω est un réel compris entre 1 et 2 [11].
Parmi les méthodes de résolution, les plus utilisées sont celles de Pour réduire le nombre d’itérations de source, on applique la
Jacobi ou de Gauss-Seidel (cf. chapitre Méthodes numériques de méthode de Wielandt [12] qui consiste à écrire l’équation (48) sous
base [A 1 220] et, pour plus d’information [11] [12] [13]. Elles peuvent la forme :
être appliquées point par point ou par lignes de mailles. θ 1–θ
A 0 – --------- P 0 Φ 0 = ------------ P 0 Φ 0 (59)
■ La méthode de Jacobi par points consiste à calculer le flux Φ dans k eff k eff
une maille en explicitant le système (48) sous la forme suivante :
(n) (n – 1) où θ est un réel légèrement inférieur à 1 [13].
A′ Φ
i = P˜ + R Φ
i i S∑ iœ œ (55) iœ Le problème est ainsi ramené au précédent avec les transforma-
œ
tions suivantes :
où A′i est la sous-matrice du système dans le point i. Dans un calcul θ
à deux groupes de neutrons, il s’agit d’un système de deux équa- A 0 → A 0 – --------- P 0 Φ 0
k eff
tions et deux inconnues. (60)
1–θ
■ La méthode de Jacobi par lignes consiste à expliciter le même P 0 → ------------ P 0 Φ 0
k eff
système pour l’ensemble des mailles i appartenant à une ligne.
Dans ce cas, il est possible de définir des matrices de transformation Une autre méthode d’accélération consiste à résoudre périodi-
qui, appliquées à l’intérieur du système, en rendent la matrice trian- quement le même problème sur un maillage plus grossier de façon
gulaire. La transformation est la suivante : à déterminer des facteurs correctifs à appliquer aux mailles du pro-
(n) (n)
blème de base. Cette méthode, très efficace lorsque le nombre
Φi = Mi + Qi + 1 ( i = i1 , i2 ) (56) élevé de points retarde trop la propagation de l’information du cal-
cul entre une maille et les autres (maillages fins), est dite coarse
où i1 et i2 sont les indices des mailles délimitant la ligne. Son intro- mesh rebalancing [13].
duction dans l’équation (55) conduit au système :
Mi = Bi + Ci Mi – 1 ( i = i1 , i2 ) (57)
4.4 Domaines d’application
où les matrices [C ] et [B ] sont calculées à partir de [A ] et [P˜ ] .
La matrice Mi pour la maille fictive précédant la maille i1 est
déduite à partir des conditions à la frontière. L’application progres- En général, le domaine d’application d’une chaîne de calcul est
sive de la relation (57) de la maille i1 à la maille i2 permet de déter- beaucoup plus sensible à la méthode d’évaluation des sections
miner les matrices M et Q dans chaque point. En parcourant ensuite efficaces (les données du calcul de diffusion) qu’à la méthode de
la ligne dans le sens contraire au précédent, les flux peuvent être calcul du flux. Les limites d’application de l’équation de la diffusion
déduits par application de la relation (56), après les avoir initialisés ont été énoncées au paragraphe 2.1. On parlera donc ici du
dans la maille fictive au-delà de i2 par prise en compte des domaine d’application de l’équation de la diffusion en différences
conditions à la frontière. finies (pour ce qui concerne l’équation de la diffusion en méthode
nodale, (§ 5.4)).
■ La méthode de Gauss-Seidel diffère de celle de Jacobi (par points
Dans les études de conception des REP, la méthode des différen-
ou par lignes) seulement par le fait que les flux utilisés pour calculer
ces finies est utilisée à une dimension pour analyser des phénomè-
la solution dans un point ou dans une ligne (à droite du signe =)
nes qui se développent principalement selon la direction axiale,
viennent en partie de l’itération en cours (ceux du point ou de la
comme les oscillations xénon et pour simuler les transitoires de
ligne déjà calculés) [11] [12]. La matrice A’ contient donc les termes
capacité de puissance. Elle est utilisée à deux dimensions (2D) en
Rij qui modelisent le couplage entre la/les maille(s) courante(s) et
géométrie plane (x, y ), pour des calculs où la composition axiale du
celles qui n’ont pas encore été calculées. Cette méthode converge
réacteur ne présente pas de discontinuités importantes. Toutefois,
plus rapidement que celle de Jacobi, mais elle a l’inconvénient de
la représentation des fuites de neutrons selon l’axe z nécessite tou-
créer une dépendance entre le calcul d’un point et celui des autres.
jours un calcul d’équivalence à partir d’une modélisation en trois
Elle est donc bien adaptée aux ordinateurs de la génération actuelle
dimensions (3D). L’utilisation de la modélisation à 2D est acceptable
mais moins adaptée aux ordinateurs à processeurs parallèles.
pour simuler certaines situations du fonctionnement courant et
Une itération du processus qui conduit à la solution du système pour les études de repositionnement du combustible des REP. Elle
(48) est appelée itération interne, une itération globale sur l’équa- n’est, en revanche, jamais satisfaisante pour les réacteurs à eau
tion (47) est appelée itération externe. Étant donné que ces proces- bouillante (REB) où le changement de phase liquide/vapeur doit être
sus peuvent être utilisés à l’intérieur d’autres algorithmes itératifs impérativement analysé en 3D.
(comme celui d’application des contre-réactions, de calcul de para-
mètre critique ou d’évolution du combustible), il est plus approprié

Toute reproduction sans autorisation du Centre français d’exploitation du droit de copie est strictement interdite.
B 3 070 − 14 © Techniques de l’Ingénieur, traité Génie nucléaire
____________________________________________________________________________________________ MÉTHODES DE CALCUL NEUTRONIQUE DE CŒUR

L’utilisation en 3D se généralise à tout type d’étude de conception : valeurs des flux moyens par maille et les valeurs des courants aux
elle s’applique tout particulièrement à l’insertion partielle des barres interfaces. Nous la rappelons ici pour faciliter la lecture :
de contrôle (la distribution axiale de sections efficaces présente alors

L
de fortes discontinuités). Les calculs à 3D de conception des REP sont
dS œ +  Σ a, i + ∑ Σ i
g → g′
∑
g g g
normalement exécutés avec au moins 9 (3 × 3) mailles par assem- J
 
Φ 0i V i
blage combustible dans les directions (x, y ) et une cinquantaine de œ=1 Sœ g′
mailles selon l’axe z, de façon à modéliser les grilles de maintien du
∑
g
combustible. – Σ g ′ → g + λ0 χ ν
g
Σ f, i
g
Φ 0i V i = 0
g′
 i 
Le passage du modèle 3D au modèle 1D se fait par des calculs
de condensation. En 2D, le niveau de détail peut être passé ■ Réévaluation des couplages
jusqu’au crayon combustible et à la lame de réfrigérant qui sépare On s’attache à définir une relation entre les niveaux de flux des
un assemblage d’un autre. En 3D, un calcul discret au niveau du mailles et les courants neutroniques. Ce sont ces relations qui cou-
crayon est toujours possible, mais il devient très volumineux. plent les nodes entre eux.
Dans les études des RNR, on adopte en général une méthodolo-
gie 3D en géométrie hexagonale z, seules les études paramétriques ■ Exemple
et de dégrossissage étant effectuées à deux dimensions en géomé- L’équation de couplage est bien souvent le résultat d’une hypo-
trie (R, z ). Le nombre de mailles utilisées pour discrétiser spatiale- thèse de continuité du flux homogène associée à une approxima-
ment le problème (aussi bien que le nombre de groupes d’énergie) tion de la loi de Fick notée Fick-gross :
conditionne la précision du calcul, par le biais de l’effet de maillage.
g g
g Φi – Φi œ
Il intervient donc directement dans la définition des marges de g
conception et le dimensionnement des provisions qui sont appli- J iœ = – D i ---------------------- (61)
hi ⁄ 2
quées aux résultats de calculs (le nombre de mailles a par ailleurs
une forte influence sur les temps de calculs). où, pour chaque couple de mailles adjacentes i et œ , on a :
g
J iœ valeur algébrique du courant neutronique dans le groupe g
passant de i à œ ;
5. Méthodes nodales g
D i valeur du coefficient de diffusion dans le groupe g
sur la maille i ;
avancées g
Φ i valeur moyenne du flux neutronique dans le groupe g
sur la maille i ;
g
Φ iœ valeur du flux surfacique dans le groupe g
Historiquement on peut dire que les méthodes nodales de à l’interface entre les mailles i et œ ;
conception sont l’aboutissement des efforts des neutroniciens pour hi largeur de la maille i.
résoudre l’équation de la diffusion. Les premières générations de
codes neutroniques souffraient des faiblesses des ordinateurs qui, Le couplage entre les mailles est assuré par les relations de
limités en taille mémoire, ne permettaient qu’une faible discrétisa- continuité des flux surfaciques et des courants :
tion spatiale. Les méthodes nodales diffusives sont le résultat de la g g g g
mise au point d’une procédure interne au processus de J iœ = J œi et Φ iœ = Φ œi (62)
convergence du code qui permet une meilleure évaluation du cou-
plage entre les mailles homogènes de calculs. Les relations (14), (61) et (62) forment le système complet
d’équations que l’on doit résoudre pour trouver la solution du pro-
Les paragraphes suivants décrivent succinctement la théorie à la blème, selon le schéma suivant :
base de cette amélioration ainsi que sa mise en place dans les
codes. Système : bilan neutronique du node n
opérateur courant : loi de Fick-gross
continuité des courants
continuité du flux homogène
5.1 Principe de l’équivalence nodale Paramètres : sections efficaces moyennes
Inconnues : flux moyens
Le problème de l’équivalence nodale a été formalisé comme suit : courants aux interfaces
supposant connue une référence neutronique hétérogène (une Le problème est donc : Quelle transformation faut-il faire subir à
solution de l’équation de Boltzmann) et supposant défini par ce système pour qu’il conduise à une solution représentative de la
ailleurs un maillage homogène, quel traitement faut-il faire subir à référence ?
l’équation de la diffusion pour que la solution trouvée sur le modèle Il est possible de reproduire sur un maillage homogène et avec
homogène soit cohérente avec la référence hétérogène ? [44]. un opérateur de diffusion les flux moyens hétérogènes, les courants
La notion de cohérence est ici fondamentale et reflète les gran- hétérogènes aux interfaces, les flux hétérogènes aux interfaces, et
deurs que l’on souhaite conserver en passant de la référence au l’ensemble des taux de réactions hétérogènes [44].
modèle homogène. Le processus d’équivalence repose sur les modifications suivan-
L’algorithme de résolution de l’équation de la diffusion est sup- tes.
posé se décomposer en deux étapes principales dont l’enchaîne- ● Sur le système
ment est le suivant.
1 ) Continuité des courants (méthode d’intégration transverse) :
■ Calcul des flux comme nous l’avons vu précédemment, l’évaluation classique des
courants d’interface se fait très souvent par l’application de la loi
On résout l’équation du bilan neutronique multigroupe sous la
de Fick en supposant que le flux varie linéairement entre la surface
forme de l’équation (14), ce qui permet d’établir un lien entre les
et le centre de la maille et que les flux surfaciques soient continus ;
cela assure une relation entre les flux moyens de deux mailles
adjacentes et le courant qu’elles échangeant.
2 ) Continuité du flux homogène : si on souhaite conserver le cou-
rant à l’interface de deux mailles adjacentes et la forme de l’opé-

Toute reproduction sans autorisation du Centre français d’exploitation du droit de copie est strictement interdite.
© Techniques de l’Ingénieur, traité Génie nucléaire B 3 070 − 15
MÉTHODES DE CALCUL NEUTRONIQUE DE CŒUR _____________________________________________________________________________________________

rateur de calcul des courants (Fick-gross ), alors il est inévitable qu’il composantes hétérogènes des facteurs de discontinuité) par les for-
faille relâcher la contrainte de continuité du flux homogène mules suivantes [46] [47] :
(figure 4). Le lien entre le flux homogène et le flux hétérogène est
introduit sous la forme de facteurs de discontinuité. À la place de la
condition de continuité des flux homogènes, on a donc la condition r, g
 motif
r, g g
Σ hét ( r ) Φ hét ( r ) d r
2

Σ hom

de continuité des flux hétérogènes. Ces facteurs dépendent de la = ------------------------------------------------------------------------
g 2
méthode d’évaluation des courants et des sections efficaces ; ils Φ hét ( r ) d r
permettent de prendre en compte l’homogénéisation de la région et motif

Φ
(63)
de compenser les approximations numériques de la méthode
g
het ( r )
d’intégration. En fait, bien souvent on distingue deux composantes dS
g S

Φ
du facteur de discontinuité : df hom = -----------------------------------------
-
g
— la composante d’hétérogénéité de la maille due à l’homogénéi- hom ( r ) dS
sation du milieu [45] ; S
— la composante numérique d’interface, qui contraint un modèle avec
r, g
Σ hom section efficace homogène pour la réaction r dans
numérique à reproduire une référence. Cette dernière est évaluée le groupe g,
à l’intérieur même du moteur nodal. r, g
Σ hét ( r ) section efficace hétérogène pour la réaction r
dans le groupe g,
g
Φ hét ( r ) flux hétérogène dans le groupe g,
g
Φ hom ( r ) flux homogène dans le groupe g,
g
df hom facteur de discontinuité pour la surface S et dans
le groupe g.

5.2 Méthode d’expansion nodale

Dans la méthode d’expansion nodale, l’opérateur de calcul du cou-

rant d’interface est basé sur la résolution de l’équation de la diffusion
selon la direction transverse à la composante du courant qu’on cher-
che.
g
Soit Φ i (u, v, w ) (u, v et w étant les coordonnées réduites) le flux
dans la maille n, on peut définir les flux transverses dans les 3
directions (u, v et w ) selon la formule suivante :


Figure 4 – Raccordement des grandeurs flux et courant
aux interfaces g g
ϕi ( u ) = Φ i ( u, v, w ) dv dw (64)
– 1 ⁄ 2 BvB1 ⁄ 2
– 1 ⁄ 2 BwB1 ⁄ 2
● Sur les grandeurs La spécificité de la méthode d’expansion nodale est d’éviter
— utilisation des sections efficaces pondérées par les flux de l’intégration analytique coûteuse des équations transverses en
référence et les volumes ; supposant que les flux 1D peuvent décomposer sur une base de
— introduction à côté des sections efficaces habituelles de fac- fonctions Fn (u ).
teurs de discontinuité (composante hétérogène). Dans ce qui suit, nous ferons l’hypothèse que la base de
Cette démarche, bien que prouvant la possibilité de réaliser un décomposition est la suivante :
algorithme de calcul homogène capable de reproduire une réfé-
rence hétérogène, ne peut en aucun cas conduire directement à une F0 ( u ) = 1
méthode opérationnelle. En effet, la construction du système F1 ( u ) = u
d’équations équivalent nécessite la connaissance de la référence
hétérogène qui n’est évidemment pas connue. Il faut donc recenser 2 1
F 2 ( u ) = u – ------
les paramètres utilisés dans le système équivalent, qui sont calculés 12
à partir de la référence et déterminer une méthode simplifiée d’éva- (65)
F 3 ( u ) =  u – --- u
2 1
luation. En pratique, on cherche la solution sur des motifs bidimen-
 4
sionnels correspondant aux zones d’homogénéisation en leur
imposant des conditions aux limites arbitraires (courant nul dans
F 4 ( u ) =  u – ---  u – ------
2 1 2 1
beaucoup de cas). On obtient les sections efficaces (y compris les  4  20

On peut donc écrire :

n=4

∑
g n
ϕi ( u ) = cg Fn ( u ) (66)
n=0

La détermination des coefficients d’expansion se fait grâce à la

résolution d’un système d’équations qui exprime que le flux est une
fonction qui doit vérifier l’équation de la diffusion, avec ses
conditions de continuité (ou discontinuité), et être en accord avec

Toute reproduction sans autorisation du Centre français d’exploitation du droit de copie est strictement interdite.
B 3 070 − 16 © Techniques de l’Ingénieur, traité Génie nucléaire
____________________________________________________________________________________________ MÉTHODES DE CALCUL NEUTRONIQUE DE CŒUR

les flux moyens par maille. Une procédure standard de détermina- ■ Flux intranodaux
tion est basée sur la résolution d’un problème dit à 2 nodes [45]. Supposons que la distribution du flux dans le groupe g puisse se
Dans cette méthode, on construit pour chaque interface entre les décomposer sur une base fonctionnelle, alors la détermination des
mailles un système d’équations qui comporte 10 G équations, G distributions se fait par la résolution d’un système d’équations. Les
étant le nombre de groupes d’énergie considéré : inconnues du système doivent être telles que les valeurs moyennes
Pour chaque groupe g, on obtient les équations suivantes : nodales soient conservées.
1) Conservation du flux moyen dans les mailles L’analyse des formes du flux homogène obtenues grâce à des
adjacentes du problème à deux nodes ....... 2 équations calculs de référence à maillage de calcul très fin conduit à la
2) Continuité du courant.................................... 1 équation conclusion qu’il est nécessaire d’avoir un développement des flux
3) Continuité du flux hétérogène...................... 1 équation à deux groupes sur une base comportant 13 fonctions. Une bonne
4) Résidus pondérés d’ordre 0, 1 et 2 pour base de représentation des flux homogènes à deux groupes est
les mailles adjacentes du problème à deux résumée ci-après ; les éléments de base sont des produits de fonc-
nodes .......................................................................____________
6 équations tions et l’on a grisé ceux qui sont retenus pour les développements
des fonctions flux.
d’où 10 équations (0)
L’ordre k de la méthode des résidus pondérés est défini par :
Groupe 1

1⁄2 2 g
g ∂ ϕi g g g g
Wk ( u ) –Di - + Σ t, i ϕ i ( u ) – S ( u ) + L ( u ) du = 0
----------- (67) P 0 (u ) P 1 (u ) P 2 (u ) P 3 (u ) P 4 (u )
2
–1 ⁄ 2 ∂u P 0 (v )
où Sg est la source neutronique et Lg les fuites transverses à la P 1 (v )
direction u.
P 2 (v )
Les fonctions de pondération W peuvent être arbitraires. En
général, on fait les choix suivants pour les trois premiers ordres : P 3 (v )
P 4 (v )
W0 (u ) = F0 (u ) W1 (u ) = F1 (u ) W2 (u ) = F2 (u )

(0)
5.3 Prise en compte de l’épuisement
Groupe 2
Les valeurs des sections efficaces utilisées dans le bilan neutro- Q 0 (u ) Q 1 (u ) Q 2 (u ) Q 3 (u ) Q 4 (u )
nique ou dans les équations de couplage 1D sont issues de calculs
de motifs en transport réalisés en amont du calcul de cœur (géné- Q 0 (v )
ration de tables paramétrées). Ces calculs représentent au mieux Q 1 (v )
les états dans lesquels les mailles de calculs se trouvent réellement
Q 2 (v )
dans le cœur au cours des cycles. La prise en compte de l’évolution
du combustible est réalisée en épuisant le motif à la puissance Q 3 (v )
moyenne du cœur.
Q 4 (v )
Dans le cas où le motif a été calculé en milieu infini, il est souvent
nécessaire d’introduire des corrections aux sections efficaces pour
prendre en compte la forme macroscopique de l’épuisement [48].
Cette modélisation est essentielle si l’on veut limiter le nombre de Les fonctions Pi sont en général choisies sous forme polyno-
mailles de calculs. miale, tandis que les fonctions Qi sont un mélange de polynômes
et de fonctions hyperboliques, ce qui permet de représenter plus
précisément les variations rapides du flux thermique aux interfaces
des assemblages.
5.4 Reconstruction fine Le calcul nodal fournit les valeurs moyennes par maille
suivantes :
(0)
Lorsqu’on effectue un calcul nodal, le processus itératif a oublié
les hétérogénéités des assemblages et il ne connaît plus que les
grandeurs moyennes par maille. Pour fournir au concepteur les Grandeurs nodales Nombre
informations crayon par crayon dont il a besoin pour ses études, il
est nécessaire de reconstruire ces informations à partir des valeurs G
g
moyennes du calcul de cœur et des calculs de motifs [49]. À titre Φm moyenne volumique
d’exemple, nous détaillerons ici la méthode de reconstruction 4G
radiale de la puissance crayon par crayon. g
Φ mS moyenne surfacique
La détermination de la puissance intranodale homogène se
déroule en plusieurs étapes qui sont : g 4G
J mS
— étape 1 : détermination des flux homogènes intranodaux ; moyenne surfacique
— étape 2 : détermination des distributions intranodales des Total 9G
sections efficaces de production d’énergie ;
— étape 3 : calcul de la distribution de puissance homogène
intranodale.
On constate que le calcul nodal ne fournit que 9 grandeurs inté-
grales. On aboutit donc à un système d’équations sous contraintes.

Toute reproduction sans autorisation du Centre français d’exploitation du droit de copie est strictement interdite.
© Techniques de l’Ingénieur, traité Génie nucléaire B 3 070 − 17
MÉTHODES DE CALCUL NEUTRONIQUE DE CŒUR _____________________________________________________________________________________________

Il faut donc définir 4 grandeurs homogènes supplémentaires pour — le calcul de cœur dispose de données sufissantes pour choisir
pouvoir résoudre le problème. dans la base de résultats du calcul de spectre les distributions de
Une possibilité, qui a fait largement ses preuves, est la détermi- puissance correspondantes (modèle de contre-réaction) ;
nation des flux homogènes aux coins des assemblages. Cette — alors, la distribution de puissance reconstruite crayon par
méthode possède l’avantage de contraindre la solution aux coins crayon s’obtient immédiatement par la formule :
des assemblages, ce qui la stabilise. La recherche des flux homo- motif
gènes aux coins se fait en estimant le flux hétérogène local. Pour Phét (u, v ) = λ Phom (u, v ) P hét ( u, v ) (69)
ce faire, on a recours à une première valeur du flux homogène au
coin. Cette procédure fournit autant d’estimations de flux aux coins où λ est une constante de normalisation définie par :


qu’il y a de mailles le partageant. Comme le flux surfacique moyen,
le flux calculé aux coins n’a aucune raison d’être continu. Seul le
flux hétérogène doit l’être ; c’est pourquoi on est amené à calculer assemblage hom
P ( u, v ) dudv
λ = --------------------------------------------------------------------------------------------------------------------

- (70)
ce dernier.
motif
Pour des raisons évidentes, et bien que le flux hétérogène soit P
assemblage hom
( u, v ) P hét ( u, v ) du dv
continu, il n’y a pas de raison pour que l’on trouve les mêmes motif
valeurs du flux hétérogène selon les mailles que l’on prend. C’est et P hét ( u, v ) est la distribution de puissance évaluée sur le motif
pourquoi, bien souvent on effectue la moyenne des valeurs trou- avec le calcul de spectre dans des conditions voisines des mailles
vées. On évalue ensuite les différentes valeurs du flux homogène du calcul de cœur.
au coin. À la fin de cette étape, nous disposons de 13 équations et
13 inconnues, ce qui permet de calculer les flux intranodaux homo-
gènes.
5.5 Domaines d’application
■ Sections efficaces intranodales
La logique de détermination est ici rigoureusement la même que
pour les flux : on suppose connue la forme du développement et, Les codes nodaux ont été développés pour traiter les cœurs de
en utilisant au maximum les résultats du calcul homogène nodal, types REP, REB ou RNR. Dans tous ces domaines ils ont permis
on se ramène à la résolution d’un système d’équations. Une bonne d’atteindre des précisions jusque-là jamais obtenues dans un
représentation des sections efficaces est donnée par une décom- contexte industriel et cela pour des gestions du combustible de
position de la section sur une base polynominale de degré 2 : plus en plus performantes et variées (gadolinium, MOX zoné, etc.).
(0) Ils permettent aussi de traiter les géométries à pas carré ou hexa-
gonal. Leur réussite est à mettre au compte de l’apport théorique
important [44] [48] qui a permis de résoudre un à un les problèmes
P0 (u ) P1 (u ) P2 (u ) rencontrés pour jeter les bases des premières chaînes nodales de
P0 (v ) conception :
— introduction d’une technique rigoureuse d’homogénéisation
P1 (v ) des constantes de diffusion (notion d’équivalence nodale) ;
P2 (v ) — introduction d’une méthode de reconstruction crayon par
crayon (déshomogénéisation ) ;
— prise en compte de l’épuisement dans les mailles de calculs
(dépendance spatiale des sections efficaces).
et le code nodal fournit les grandeurs suivantes :
(0) Depuis, de nombreux codes de calculs basés sur ces principes
de base ont été développés et on assiste à la floraison de nom-
breuses chaînes de calculs presque identiques pour ce qui est des
Grandeurs nodales Nombre principes de base. Ainsi, tous les grands concepteurs ont choisi de
réaliser leurs chaînes nodales sur l’enchaînement d’un code de cel-
≡≡≡≡≡≡ lule d’assemblage (le motif d’homogénéisation est un assemblage)
g
κΣ f ,m moyenne volumique G
et d’un code nodal avancé permettant de reconstruire entre autres
======= la puissance crayon par crayon.
g
κΣ f ,mS moyenne surfacique 4G La seule chose qui différencie encore les codes nodaux, c’est le
lien entre le calcul de cœur et les calculs d’assemblages, c’est-à-dire
Total 5G la modélisation des contre-réactions. Ces contre-réactions permet-
tent au calcul de déterminer les grandeurs d’assemblages dont il a
besoin (sections efficaces, distribution de puissance, etc.) pour
Contrairement au cas des flux, nous disposons directement du représenter au mieux les grandeurs correspondantes pour chaque
nombre correct de données pour reconstruire la fonction intrano- maille de la géométrie du cœur.
dale correspondant à la section de production d’énergie. Quel que soit le modèle, tous essaient de prendre au mieux les
effets d’historique de l’épuisement et d’environnement qui ne peu-
■ Puissance homogène intranodale vent être vus par des calculs en milieu infini. En effet, le dévelop-
La puissance homogène est obtenue directement par le produit pement d’assemblages de plus en plus particuliers avec notamment
des distributions de flux et de sections : des hétérogénéités aux frontières (gadolinium périphérique, zonage
MOX, etc.) et les gestions du combustible avec des grappes forte-
G
ment insérées dans le cœur rendent nécessaire l’introduction de
∑
g g
P hom ( u, v ) = κ Σ f ( u, v ) Φ ( u, v ) (68) raffinements à la méthode [50].
g=1
Les méthodes nodales ont aussi été appliquées avec succès en
Si l’on suppose que : cinétique pour traiter des transitoires accidentels, ce qui fait que
— le calcul de cellule fournit, en plus des sections efficaces, des les concepteurs ont à leur disposition toute une panoplie d’outils
distributions de puissance crayon par crayon sur son motif de permettant de mener à bien l’ensemble des études de neutronique
calcul ; dont ils ont la charge. Les méthodes nodales sont maintenant arri-

Toute reproduction sans autorisation du Centre français d’exploitation du droit de copie est strictement interdite.
B 3 070 − 18 © Techniques de l’Ingénieur, traité Génie nucléaire
____________________________________________________________________________________________ MÉTHODES DE CALCUL NEUTRONIQUE DE CŒUR

vées à maturité mais des développements restent à envisager pour sa contribution asymptotique à la réaction en chaîne. Ainsi, à un
le futur, notamment en ce qui concerne : neutron de très haute énergie situé près du bord du cœur correspond
— le raffinement du traitement des équations transverses ; une faible valeur du flux adjoint, alors qu’un neutron thermique au
— l’optimisation de l’algorithmique des codes ; centre de la région combustible est caractérisé par une valeur élevée
— l’amélioration du traitement des conditions aux frontières de cette fonction.
(réflecteur). L’existence du flux adjoint trouve sa justification dans la direction
La prochaine grande étape industrielle serait l’abandon de l’équa- privilégiée qui conditionne le déroulement des histoires neutro-
tion de la diffusion et le développement d’outils capables de résou- niques en physique des réacteurs. Les neutrons possèdent à la nais-
dre directement l’équation du transport en 3D. sance une énergie élevée (2 MeV, en moyenne, dans les REP
classiques, à peine plus dans les RNR) mais, lorsqu’ils sont absorbés
ou quittent le système, ils n’en conservent qu’une partie réduite.
Pendant leur vie à l’intérieur du cœur, ils cèdent, en effet, tout ou
6. Méthodes de perturbation partie de leur énergie cinétique à l’environnement. Ceux qui pro-
duisent des fissions ont des descendants dont l’énergie est plus éle-
en physique des réacteurs vée que la leur, et ainsi de suite. Il s’instaure de cette façon un
processus cyclique en énergie qui se traduit en forme mathématique
par la non-réversibilité de certains opérateurs figurant dans l’équa-
tion d’état du système.
6.1 Méthode CPT (Classical Perturbation
■ Usachev en 1955 [52] proposa la première formulation complète
Theory ) et cohérente de la CPT, en y incluant aussi le calcul du temps de vie
des neutrons prompts et la fraction efficace de neutrons retardés ;
Les méthodes de perturbation ont joué un rôle important en phy- (pour plus d’information sur ces paramètres essentiels en cinétique
sique des réacteurs depuis les débuts des applications de l’énergie des réacteurs, cf. chapitre Théorie des réacteurs nucléaires [B 3 025]
nucléaire, notamment dans le domaine des expériences critiques et [2] [7]).
et des mesures. En accord avec cette formulation, toute variation d’une
composante générique des opérateurs de l’équation de Boltzmann
■ Déjà Wigner en 1945 [51] en proposait la formulation plus élé- du système produit une variation de la réactivité qui vaut * :
mentaire, la CPT (Classical Perturbation Theory, théorie classique
des perturbations), issue directement des concepts classiques de la *
〈 Φ 0 , [ δA + ( 1 – ρ ) δP ] Φ ′〉
mécanique quantique, pour mesurer le poids en réactivité de petits δ ρ = ---------------------------------------------------------------------
*
- (74)
échantillons de matériaux introduits dans une pile atomique. 〈 Φ 0 , P ′ Φ ′〉
Dans les premières années du développement du nucléaire civil,
les méthodes de perturbation occupèrent une place très impor- avec δρ variation de réactivité recherchée,
tante dans les analyses de sûreté car, en raison des moyens de cal- Φ’ solution de l’équation d’état du système dans les
cul très limités, les études étaient généralement effectuées dans conditions asymptotiques après perturbation,
l’approximation à zéro dimension, moyennant l’utilisation de coef- [δA + (1 – ρ) δP ]
ficients de sensibilité obtenus justement à l’aide de la méthode représentant la perturbation générique des opérateurs
CPT. de l’équation de Boltzmann.
La définition du flux adjoint, parfois appelé aussi importance neu- (*) Cette équation est obtenue, sans introduire aucune approximation, à partir de
tronique s’est avérée essentielle au développement de la méthode l’équation d’état du système après perturbation qui s’écrit de façon générique :
CPT. A’ Φ’ + λ’ P’ Φ’ = 0
Le flux adjoint est la solution unique de l’équation (11) : avec A’ = A0 + δA,
λ’ = 1 ⁄ k eff
′ = (1 – ρ’) avec ρ’ = ρ + δρ,
t * t *
A0 Φ0 + λ0 P0 Φ0 = 0 (71) P’ = P0 + δP,
Φ’ = Φ0 + δΦ
assortie de conditions aux frontières génériques du type : A0 , P0 , Φ0 , ρ étant respectivement les opérateurs, la fonction propre et la réactivité de
l’état de référence avant perturbation.
*
* ∂Φ 0 En remplaçant directement les définitions dans l’équation, on a :
α Φ 0 + β ---------- = 0 (72)
∂n (A0 + δA) Φ’ + [1 – (ρ + δρ )] (P0 + δP ) Φ’ = 0

avec α, β constantes de proportionnalité et n direction générique. En multipliant à gauche par Φ *0 , on obtient :

Cette équation est l’adjointe de celle du bilan neutronique du < Φ *0 , ( A 0 + δA ) Φ ′> + < Φ *0 , [ 1 – ( ρ + δ ρ ) ] ( P 0 + δP ) Φ ′> = 0
système en forme homogène [équation (10)] :
et :
A 0 Φ 0 + λ0 P0 Φ 0 = 0 < Φ *0 , A 0 Φ ′> + < Φ *0 , ( 1 – ρ )P 0 Φ ′> + < Φ *0 , δA Φ ′>
+< Φ *0 , ( 1 – ρ ) δP Φ ′> – < Φ *0 , δ ρ P Φ ′> = 0
Les flux adjoint possède des propriétés mathématiques similai-
res à celles du flux. Il est notamment associé à la même valeur pro- Réarrangeant et faisant appel aux propriétés du flux adjoint, on obtient :
pre λ0 et défini positif sur tout le domaine énergétique et spatial.
De plus, il satisfait les conditions d’orthogonalité suivantes : < Φ *0 , [ δA + ( 1 – ρ )δP ] Φ ′> = < Φ *0 , δ ρ P′ Φ ′>

ce qui donne, enfin :

〈 Φ*0 , P0 Φ0 〉 = C  〈 Φ *0 , [ δA + ( 1 – ρ ) δP ] Φ ′〉
 quel que soit i ≠ 0 (73) δ ρ = -----------------------------------------------------------------------
〈 Φ *0 , P ′Φ ′〉
-
〈 Φ*
0 , P0 ψ i 〉 = 0 

ψi étant les harmoniques de l’opérateur de Bolzmann, ainsi,

ψ0 = Φ 0 .
Le flux adjoint possède une signification physique originale car
il mesure l’importance de chaque neutron du système vis-à-vis de

Toute reproduction sans autorisation du Centre français d’exploitation du droit de copie est strictement interdite.
© Techniques de l’Ingénieur, traité Génie nucléaire B 3 070 − 19
MÉTHODES DE CALCUL NEUTRONIQUE DE CŒUR _____________________________________________________________________________________________

Si on fait l’approximation que la perturbation n’affecte pas (ou qui respecte strictement la condition d’orthogonalité :
peu) le flux, on peut simplifier l’équation précédente qui devient
l’expression perturbative classique au premier ordre : 〈 ∗ , Φ 0〉 = 0 (79)

〈 Φ* 0 , [ δA + ( 1 – ρ ) δP ] Φ 0〉 Cela permet d’écrire l’équation (77) sous la forme :

δ ρ = -------------------------------------------------------------------------
- (75)
〈 Φ* 0 , P 0 Φ 0〉
δ 〈 δh 1 , Φ 0〉 〈 δh 2 , Φ 0〉
- – -------------------------- + 〈 ∗ , δ Φ〉
-------- = ------------------------- (80)
Les équations (74) et (75) peuvent être simplifiées ultérieurement  〈 h 1 , Φ 0〉 〈 h 2 , Φ 0〉
sur la base d’hypothèses variées, qui vont de l’approximation zéro
dimension (milieu homogène infini) à celle monocinétique (un seul Il s’agit maintenant de rechercher une solution de l’équation
groupe d’énergie), selon le type de problème à traiter et la nature adjointe avec ∗ comme terme inhomogène, à savoir :
des études. t t
A 0 ψ ∗ + λ 0 P 0 ψ ∗ = ∗ (81)
On peut, notamment, en profitant de la propriété transitive, utiliser
des formulations dérivées de celles-ci pour évaluer la contribution
En accord avec l’alternative de Fredholm (cf. chapitre Équations
individuelle de chaque isotope à l’écart en réactivité entre deux états
aux dérivées partielles [A 650] dans le traité Sciences fondamentales
d’un système, comme il est fait par le sous-programme Chloé du
et [1], cette équation admet une solution unique, connue avec un
code Apollo-II [53].
degré de liberté, sous la forme :

ψ∗ = ψ + a Φ 0 (82)
6.2 Méthode GPT (Generalized Afin de garantir l’unicité de la solution, il faut choisir le coefficient
Perturbation Theory ) a de façon à faire disparaître le deuxième terme de l’équation. Cette
opération de filtrage, effectuée systématiquement pendant le calcul
de ψ*, est connue sous le nom de décontamination du mode fon-
Au milieu des années 1960, Usachev [54], Lewins [55] et damental.
Gandini [56] généralisèrent la théorie des perturbations à toute
sorte d’observable, paramètre physique pouvant être mesuré (et/ou En remplaçant maintenant l’équation (81) dans la partie spectrale
calculé), qui peut être représenté sous la forme d’un rapport de de l’équation (77), on obtient, au premier ordre :
quantités intégrales, linéaires du flux et/ou du flux adjoint.
δ spectral
= – 〈 ψ∗ , ( δA + λ 0 δP ) Φ 0〉 + δ λ 0 〈 ψ∗ , A 0 Φ 0〉 (83)
Pour reconstruire la valeur d’un observable d’un système multi- --------

plicateur dans des conditions perturbées, la théorie généralisée des
perturbations, ou GPT (Generalized Perturbation Theory ), appelée Si, conformément aux hypothèses, la perturbation ne modifie pas
aussi, historiquement, théorie des perturbations généralisée, fait la réactivité du système, le terme δ λ0 <ψ *, A 0 Φ0> est nul ; le cas
appel à une fonction particulière, dite fonction importance généra- échéant, il peut être pris en compte comme correction dite de reactivity
lisée, qui est la solution unique de l’équation du bilan adjointe du reset. Cette correction (mise à zéro de la réactivité = rétablissement
système à laquelle on ajoute une terme inhomogène, appelé source des conditions de criticité du système après perturbation) simule
généralisée. Ce terme, qui est strictement orthogonal au mode l’action des moyens de contrôle utilisés dans le but de rétablir la
fondamental, contient, au premier ordre, l’information sur la nature criticité du système, après perturbation. Il existe plusieurs façons
de la perturbation et son impact sur le flux. d’introduire cette correction selon le mode opératoire que l’on veut
Soit  un observable générique écrit sous la forme d’un rapport simuler [57] [58].
de quantités intégrales linéaires dans le flux :
〈 h 1 , Φ 0〉
 = ----------------------
- (76) 6.3 Développements récents
〈 h 2 , Φ 0〉

où h1 et h2 sont deux fonctions positives, par exemple des sections

Les années 1970 et le début des années 1980 ont marqué un
efficaces.
nouvel essor de la théorie GPT avec, notamment :
Si l’on introduit dans le système une perturbation qui ne modifie — son extension au domaine non linéaire [59] ;
pas, ou peu, la réactivité, la variation relative de la fonctionnelle (76) — sa première application aux problèmes d’évolution du
s’exprime par : combustible [60] ;
— les travaux de Cacuci visant le développement et la systéma-
δ 〈 δh 1 , Φ 0〉 〈 δh 2 , Φ 0〉 〈 h 1 , δ Φ〉 〈 h 2 , δ Φ〉 tisation de la formulation différentielle sur une assise mathématique
-------- = -------------------------- – ---------------------------- + ------------------------- – -------------------------
 〈 h 1 , Φ 0〉 〈 h 2 , Φ 0〉 〈 h 1 , Φ 0〉 〈 h 2 , Φ 0〉 rigoureuse [61] et ceux de Gandini, qui ont abouti à la méthode dite
(77) heuristique, basée sur le principe de la conservation de l’importance
neutronique [62].
↑ ↑
Un autre volet des méthodes de perturbations a été consacré à
effet direct effet spectral l’étude des variations d’interface géométrique des milieux pour
évaluer, par exemple, les conséquences sur la réactivité des mou-
Le calcul de l’effet direct est effectué à l’aide de la méthode CPT,
vements des grappes de contrôle, des fluides ou du combustible.
alors que celui de l’effet spectral peut être fait, au premier ordre en
δΦ, à l’aide de la méthode GPT. Aux États-Unis, Pomranging [63] et Ranhema [64] ont proposé une
méthode basée sur le développement explicite de la valeur propre
On définit alors une source généralisée ∗ sous la forme :
en deux termes, l’un géométrique et l’autre lié au bilan matière. En
h1 h2 France, Porta [65] a mis au point une méthode axée sur le dévelop-
∗ = ----------------------
- – ----------------------
- (78) pement implicite de la valeur propre du système conditionné, où le
〈 h 1 , Φ 0〉 〈 h 2 , Φ 0〉 conditionnement numérique du problème le ramène à une résolu-
tion de type GPT.
Récemment, l’attention des chercheurs s’est portée sur un groupe
de méthodes qui recherchent plus particulièrement la reconstruc-

Toute reproduction sans autorisation du Centre français d’exploitation du droit de copie est strictement interdite.
B 3 070 − 20 © Techniques de l’Ingénieur, traité Génie nucléaire
____________________________________________________________________________________________ MÉTHODES DE CALCUL NEUTRONIQUE DE CŒUR

tion explicite de la forme du flux dans les conditions asymptotiques

après perturbation. Gandini, le premier, a proposé ces méthodes en 7. Schémas de calcul
1978 [66] en préconisant un développement limité du flux perturbé
sur la base des fonctions propres de l’opérateur de Boltzmann du
système de référence, non perturbé. Cette méthode, dite Méthode 7.1 Structure et qualification
Standard (SM), s’est avérée performante dans des cas simples, non
symétriques et à un seul groupe d’énergie, mais elle n’a jamais
des chaînes de calcul
donné pleine satisfaction dans des cas plus généraux, notamment
multigroupes. Cela s’explique par le fait que le jeu de fonctions pro- La conception neutronique d’une centrale nucléaire de puissance
pres de type λ, retenues comme base pour la reconstruction, n’est comporte l’exécution d’un grand nombre de calculs dont le but est
pas complet par rapport à la distribution énergétique des neutrons*. d’appréhender les grandeurs qui décrivent son comportement,
Pour résoudre cette difficulté, Gandini a proposé l’utilisation du jeu dans toutes les conditions de fonctionnement. Nous avons vu que,
(complet, cette fois) des fonctions propres de type ω, issues de en raison de la complexité des équations et du nombre de para-
l’opérateur inhomogène de la cinétique [67]. Cette approche s’est mètres en jeu, ces calculs ne peuvent pas être effectués par voie
avérée d’application difficile, en raison du fait que toutes les valeurs analytique : les études de conception nécessitent donc des codes
propres de type ω sont complexes, à l’exception de la première. capables de résoudre les problèmes par voie numérique. Ces
(*) L’hypothèse selon laquelle la distribution énergétique des neutrons issus de fission codes sont utilisés sur des ordinateurs de capacité élevée, que ce
est indépendante des énergies de ceux qui engendrent les fissions est à l’origine de
l’incomplétude du jeu des fonctions propres de l’opérateur de Boltzmann, dans les cas
soient des machines de grande taille ou des stations de travail
multigroupes. Il suffit d’imposer un très faible couplage entre ces deux familles de neu- connectées en réseau.
trons, comme suggéré dans [71] pour rendre complet le je de ces fonctions propres.
Les codes, à l’image de l’approche retenue dans la conception
Pour pallier les inconvénients et les limitations de la méthode SM, industrielle des tranches nucléaires, sont conçus dans le but précis
de nombreuses méthodes ont été proposées, notamment dans le de résoudre, chacun de façon indépendante des autres, une partie
but de compléter le jeu de fonctions propres de l’opérateur de Boltz- spécifique et limitée des problèmes à traiter et de fournir, en sortie,
mann**. On retiendra la méthode des pseudo-harmoniques, déve- des résultats utilisables comme données d’entrées pour les codes
loppée par EDF dans les années 1985 [68], et la méthode French, situés en aval. Il en découle que, dans leur ensemble, les différents
proposée par le CEA [69]. codes constituent les macromodules (ou process ) d’une structure
(**) Dans les systèmes symétriques, en faisant abstraction des problèmes liés à l’énergie modulaire rendue homogène et cohérente par les liens existant
déjà mentionnés, la base des fonctions propres de l’opérateur de Boltzmann peut s’avérer entre les éléments qui la composent. Souvent, cette structure
incomplète par dégénérescence. En raison de la méthode adoptée pour séparer les diverses
fonctions propres, méthode dite de puissance (la même qui est utilisée couramment dans complexe profite d’un macro-langage d’utilisation commun à tous
les calculs du flux) assortie d’une filtrage systématique des modes inférieurs par ortho- les modules, d’une gestion dynamique de la mémoire et de procé-
gonalisation, la base résultante est incomplète, si le système est symétrique, en raison de dures d’entrée-sortie et gestion de l’information standardisées. On
la dégénérescence des valeurs propres des harmoniques symétriques par rotation [70].
appelle l’ensemble une chaîne de calcul de conception.
Tout récemment, de nouveaux progrès dans ce domaine ont été
Chaque code, en général, présente une panoplie d’options de cal-
enregistrés grâce à Framatome, qui a proposé une façon originale
cul correspondant aux niveaux différents d’approximation adoptés
de généraliser la méthode SM aux problèmes symétriques et aux
dans la description du comportement du réacteur. Le choix d’un
cas multigroupes [70] [71].
ensemble homogène d’hypothèses parmi celles disponibles
Ce paragraphe ne peut que présenter un très bref aperçu de la dépend alors du type de système physique que l’on veut étudier
méthodologie. Pour une information plus complète, nous renvoyons et du degré de précision souhaité dans sa représentation. Le choix
à la littérature et, tout particulièrement, aux ouvrages de compilation des options doit être effectué en respectant strictement la cohé-
sur le sujet [14] [15] [16]. rence entre les diverses étapes du calcul, de façon à éviter que la
Nombreuses sont les applications des méthodologies perturba- précision atteinte à un moment donné dans la modélisation d’un
tives aux études de conception des réacteurs. Dans le domaine des phénomène particulier ne soit remise en cause par d’éventuelles
REP, deux parmi les plus récentes méritent d’être signalées tout approximations plus grossières introduites par la suite. Cet ensem-
particulièrement, à savoir : ble cohérent et homogène d’approximations, de choix et d’options
— l’analyse de la sensibilité de la distribution de puissance aux constitue ce que l’on appelle, en général, une procédure de calcul
incertitudes de fabrication du combustible. La méthodologie pro- de conception.
posée utilise la méthode GPT pour évaluer la matrice de sensibilité Le processus de qualification d’un schéma de calcul comporte
des distributions locales de puissance aux paramètres de projet. Un l’exécution d’une grande quantité de comparaisons calcul/mesure,
couplage avec l’analyse statistique permet ainsi très facilement portant sur un nombre restreint de paramètres intégraux repré-
d’apprécier l’impact des incertitudes de fabrication sur la puissance sentatifs du comportement du réacteur et des modélisations rete-
locale du cœur. La méthode, qui a atteint le stade préindustriel, a nues du fonctionnement du réacteur sur la base d’un ensemble
fait l’objet de publications internationales, notamment [72] ; d’informations expérimentales représentatif de la multiplicité des
— l’optimisation de la gestion du combustible. Le sujet a été situations se présentant en exploitation* [87] [88].
abordé dans le cadre du développement du code Formosa, qui (*) Le processus de qualification d’une chaîne de calcul neutronique dépend étroitement
applique les méthodes de perturbation à la recherche de des options de calcul qui la caractérisent. Il comporte, en tout cas, un certain nombre d’éta-
pes distinctes qui visent chacune à qualifier un paramètre intégral de façon indépendante
configurations optimales du cœur en les couplant avec une des autres. C’est le cas, par exemple de la réactivité, et de la distribution de puissance.
technique de minimisation performante (recuit simulé). Une docu- Le remontage des incertitudes se fait en tenant compte de la démarche adoptée pendant
mentation assez détaillée sur le sujet, assortie d’applications peut la qualification. Cela explique le lien de dépendance étroite existant entre les valeurs des
incertitudes et la qualification des chaînes.
être trouvée dans [73] [74].
L’information que l’on tire de ces comparaisons permet d’attein-
dre les deux objectifs principaux de la qualification :
— d’une part, la définition précise du domaine d’applicabilité du
schéma de calcul, c’est-à-dire la détermination des systèmes qu’il
est convenable d’étudier sans modification (ce domaine est défini
en termes d’intervalles de variabilité des paramètres fondamen-
taux, comme, par exemple, dans les REP classiques, le rapport de
modération, l’enrichissement du combustible et/ou le type d’absor-
bants utilisé dans les barres de contrôle) ;

Toute reproduction sans autorisation du Centre français d’exploitation du droit de copie est strictement interdite.
© Techniques de l’Ingénieur, traité Génie nucléaire B 3 070 − 21
MÉTHODES DE CALCUL NEUTRONIQUE DE CŒUR _____________________________________________________________________________________________

— d’autre part, la quantification du niveau d’incertitude affectant Le système Smart est basé sur une modélisation 3D de type nodal
l’évaluation de tels paramètres dans le cas où la procédure prévue avancé, à deux groupes d’énergie. Il met en œuvre la méthode
a été appliquée conformément aux prescriptions, dans l’intervalle d’expansion nodale utilisant des coefficients de discontinuité intra-
de validité défini ci-dessus. nodaux et des formes intranodales polynomiales pour les fuites
L’information expérimentale utilisable dans le processus de qua- transverses, les sections efficaces, etc. Ce système est complété par
lification peut avoir son origine, selon le type de paramètres une modélisation 1D axial.
recherchés et la quantité d’informations disponibles : Le modèle de contre-réactions de Smart qui intègre un calcul de
— dans les expériences spécifiques et ciblées, effectuées dans thermohydraulique canal fermé est basé sur une tabulation multi-
les maquettes critiques des centres de recherche ; paramétrée des sections efficaces. Les paramètres utilisés sont
— dans les essais de démarrage et le retour d’exploitation du l’épuisement du combustible, la concentration en bore, le niveau
parc nucléaire en fonctionnement, qui constituent le retour d’expé- xénon, la densité de l’eau, la température combustible, la présence
rience d’exploitation. de grappes de contrôle ainsi qu’un paramètre représentatif de l’his-
toire spectrale du combustible. L’évolution du combustible est cal-
Les expériences critiques fournissent des informations très pré- culée avec un modèle d’épuisement microscopique. La chaîne
cises** sur les valeurs de la réactivité et les distributions fines des d’évolution comprend, pour une utilisation standard, 36 noyaux
taux de réaction. Elles permettent, en outre, d’appréhender les lourds et produits de fission.
effets d’interaction des absorbants et d’en mesurer l’efficacité rela-
tive. Le retour d’exploitation porte, pour l’essentiel, sur les cartes La reconstruction des valeurs locales de flux, de puissance, de
de flux prises par scrutation à des instants précis de la vie du sys- taux de réaction et d’épuisement du combustible est réalisée par la
tème, les durées de cycle, les efficacités différentielles des groupes convolution d’un flux homogène intranodal et de facteurs de forme
d’absorbants, les effets de puissance et température. Il fournit, en hétérogènes tabulés lors de calculs de réseau en transport (§ 5).
outre, les informations permettant d’effectuer, sous certaines Cette méthode de reconstruction de puissance est également utili-
conditions, des analyses de combustible irradié. sée pour évaluer l’activité au niveau des tubes instrumentés pour le
(**) La précision expérimentale obtenue sur des maquettes critiques particulièrement traitement des cartes de flux.
soignées peut être assez élevée, notamment en ce qui concerne les distributions de taux Pour rendre l’utilisation de la chaîne SCIENCE conviviale et effi-
de fission dans les systèmes hétérogènes. Citons ici, pour mémoire, les travaux effectués
dans le cadre de la mise en œuvre de la maquette Épicure dans le réacteur Éole du Centre cace, celle-ci est dotée d’une environnement d’utilisation moderne
d’études nucléaires de Cadarache [85] [86]. s’appuyant sur les techniques d’interfaces homme/machine dispo-
La qualification est indispensable à l’obtention des autorisations nibles sur station de travail. Elle se présente en fait comme un bureau
requises des Autorités de Sûreté. Elle est tout aussi nécessaire d’étude complet pour l’ingénieur et lui apporte tous les services dont
pour apprécier la possibilité de traiter des modifications de projet il a besoin [76].
ou des demandes de l’exploitant, qui conduisent à l’utilisation de Ainsi, l’utilisateur de SCIENCE à sa disposition un éditeur graphi-
la chaîne de calcul en dehors du domaine d’utilisation déjà qualifié. que qui lui permet de décrire l’étude qu’il veut réaliser sous la forme
d’un enchaînement de tâches élémentaires. Les modules qui exé-
cutent ces tâches sont choisis dans une base de procédures de cal-
cul fournie à l’utilisateur. La saisie des données est réalisée grâce
7.2 Chaîne de calcul de conception à des interfaces graphiques associées à chaque procédure. L’éditeur
des REP graphique général permet également de lancer un calcul, ou un
enchaînement de calculs, et d’en contrôler l’exécution.
Nota : se reporter utilement aux chapitres spécialisés sur les REP [B 3 000] [B 3 060] Les principaux résultats physiques peuvent être repris par un
dans ce traité. outil de posttraitement : l’utilisateur visualise alors ses résultats sur
La chaîne de calcul SCIENCE (Système de Calcul Intégré pour des écrans graphiques ou peut en extraire certains pour les trans-
l’Étude Neutronique des Chaudières à Eau), mise en service en férer à des outils classiques de traitement de données.
1995 à Framatome pour la conception neutronique des REP, est un La chaîne SCIENCE intègre également une base de données rela-
bon exemple d’un outil moderne d’analyse du cœur qui allie une tionnelle destinée à gérer les données technologiques correspon-
modélisation physique de haut niveau et un environnement d’uti- dant aux différents projets et les informations concernant les
lisation évolué [75]. études déjà réalisées.
La qualité et la précision des résultats obtenus avec cette chaîne La qualification de la chaîne SCIENCE est issue essentiellement
reposent en tout premier lieu sur la qualité des modélisations phy- de l’exploitation du retour d’expérience des réacteurs de puis-
siques qu’elle met en œuvre : sance. En effet, le suivi d’exploitation accumulé sur l’ensemble des
— des calculs d’assemblages en transport, par la méthode des réacteurs Framatome en fonctionnement ainsi que sur les autres
probabilités de collision, réalisés au moyen du code Apollo 2-F (ver- réacteurs pour lesquels Framatome a fourni des recharges repré-
sion industrielle du code Apollo-II ). Le code Apollo-II est développé sente plus de 400 cycles-réacteurs. Cet ensemble couvre une très
par le CEA ; il a été adapté, industrialisé et qualifié par Framatome large gamme de gestions du combustible : taux d’épuisement, ges-
pour ses besoins propres et il est intégré à SCIENCE sous l’appel- tions classiques par tiers, cycles longs, gestions à haut burnup,
lation Apollo-2-F ; gestions avec combustible gadolinié ou avec combustible pluto-
— des calculs de cœur en 3D par la méthodologie nodale avec nium, etc.
méthode d’expansion et reconstruction fine de puissance couplée Les essais physiques de démarrage et les cartes de flux pério-
avec un modèle de contre-réactions par tables multi-paramétrées, diques fournissent de nombreux éléments de validation des claculs
réalisés au moyen du code Smart. en termes de réactivité. À titre d’exemple, la figure 5a présente un
Le schéma de calcul mis en œuvre avec Apollo-2-F pour la généra- dépouillement, sur environ 200 cas, des mesures de concentration
tion des grandeurs nécessaires à un calcul Smart est basé sur le en bore critique réalisées lors de chaque carte de flux. Les résultats
découplage du calcul de cellule et du calcul spatial (calcul de cœur ). montrent une parfaite adéquation des calculs.
Le calcul de cellule est réalisé avec un découpage énergétique à 99 L’analyse de la qualité des distributions de puissance calculées
groupes sur une géométrie assemblage simplifiée par des regroupe- avec Smart est effectuée sur la base de comparaisons directes sur les
ments (modélisation multicellule ), chaque cellule étant discrétisée. activités dans les tubes instrumentés. La figure 5b correspondant au
Le calcul spatial est réalisé avec un découpage énergétique réduit à 6 traitement d’environ 150 cartes de flux concerne les activités inté-
groupes sur une géométrie assemblage discrétisée comportant des grées sur la hauteur active du cœur et la figure 5c concerne les acti-
cellules homogènes. vités locales calculées en 3D (environ 100 000 points). Elle permet de

Toute reproduction sans autorisation du Centre français d’exploitation du droit de copie est strictement interdite.
B 3 070 − 22 © Techniques de l’Ingénieur, traité Génie nucléaire
____________________________________________________________________________________________ MÉTHODES DE CALCUL NEUTRONIQUE DE CŒUR

— une méthodologie permettant de passer du calcul projet à la

situation réelle (inaccessible au calcul direct) par évaluation de cor-
rections de méthode ;
— des incertitudes induites par la méthodologie précédente et les
données de base.
Les données de base sont pour l’essentiel les sections efficaces
de chaque élément du cœur (nuclides lourds, constituants de
l’acier, absorbants, sodium).
Le fichier de données nucléaires est défini dans un découpage fin
à 600 groupes d’énergie couvrant la gamme depuis l’énergie ther-
mique jusqu’à 14,5 MeV. Seules les sections efficaces et les données
dont les variations avec l’énergie peuvent être considérées comme
lentes sont tabulées dans un découpage plus large à 25 groupes
d’énergie.
Les sections efficaces du formulaire Carnaval IV sont ajustées :
les valeurs des données de base sont modifiées pour minimiser les
écarts calculs-expérience relatifs à un grand nombre de paramètres
neutroniques intégraux et ponctuels. Les données expérimentales
utilisées pour cet ajustement proviennent de mesures soit sur
maquette critique, soit sur réacteur de puissance.
Les sections efficaces nécessaires pour les calculs de projet sont
générées par le code de cellule Hétaire. Ce code est un exemple de
l’optimum recherché au niveau du calcul du bilan neutronique entre
la rapidité du calcul (≈ 10 s par cellule sur Cray-YMP ) et la précision
des résultats. Le code Hétaire résout l’équation intégrale du trans-
port sur la géométrie du crayon combustible par une méthode de
probabilités de collision en géométrie hétérogène plane ou cylin-
drique, dans l’approximation B1 (le remplacement du code Hétaire
par Ecco est envisagé à plus long terme).
Les sections efficaces ainsi calculées tiennent compte, dans le
combustible, de l’hétérogénéité de l’assemblage (sections efficaces
homogènes équivalentes ).
Après cette étape, on aborde le calcul du flux macroscopique à
l’aide des codes mentionnés ci-après.
L’ensemble de programmes de calcul permettant, dans le cadre
d’un projet de RNR, de déterminer les caractéristiques neutroniques
d’un réacteur à un instant donné ainsi que l’évolution de ces carac-
téristiques en fonction de l’irradiation a été rassemblé dans le pro-
giciel modulaire CCRR (Calcul de Cœurs de Réacteurs Rapides ). La
version actuelle de CCRR (la 5.3) est aussi la version initiale d’Eranos,
le système européen commun de programmes de calcul, cœur et
radioprotection, des RNR [78] [79].
Les modules de CCRR permettent de réaliser les calculs inhérents
à toutes les étapes d’une étude neutronique de cœur.
Figure 5 – Distribution des écarts calcul/mesure L e s c a l c u l s d e b a s e ( r é a c t i v i t é d u c œ u r, e f fi c a c i t é d e s
sur la réactivité, les activités intégrées et locales 3D absorbants [80], distributions de puissance, etc.) sont effectués, en
trois dimensions et en géométrie hexagonale. L’équation de la diffu-
sion des neutrons est résolue par la méthode des différences finies.
valider le calcul du facteur d’élévation d’enthalpie pour lequel une Les calculs en géométrie hexagonale à deux dimensions sont utili-
précision meilleure que 4 % est atteinte. sés, principalement, pour déterminer certaines corrections de
On trouvera les principaux éléments de qualification de la chaîne méthode et, dans les études préliminaires, pour déterminer les flux
SCIENCE dans [77]. d’essais nécessaires pour la résolution de l’équation de la diffusion
des neutrons par la méthode de synthèse (en géométrie hexagonale
3D). Les calculs en géométrie cylindrique permettent de déterminer
des flux de condensation, des facteurs correctifs et des coefficients
7.3 Formulaire de projet des RNR de réactivité et de contre-réaction thermique. Le module de CCRR
Bistro, qui exécute des calculs 2D en géométrie (R, z ) et (x, y ) dans
l’approximation du transport, est utilisé principalement pour déter-
Nota : se reporter utilement aux chapitres spécialisés sur les RNR dans ce traité
[B 3 061] [B 3 063] [B 3 170].
miner des facteurs correctifs (diffusion-transport, homogène-hétéro-
gène) et dans les études de radioprotection.
À présent, les calculs neutroniques de projet RNR sont effectués
à l’aide du formulaire Carnaval IV développé par le CEA. Le concepteur est confronté au problème de ne pas pouvoir
représenter exactement la géométrie réelle de la configuration exa-
Ce formulaire comporte : minée. En effet, cela supposerait, au moins, un calcul en théorie du
— des données de base (sections efficaces microscopiques, don- transport, 3 dimensions, 25 groupes d’énergie, maillage infini
nées de décroissance, compositions isotopiques, stoechiométri- (espace continu) avec représentation exacte des éléments non
ques, etc.) ; fissiles (absorbants et diluants). Aucun code de calcul ne répond à
— des codes de calcul ; la question, et quand bien même il existerait, son coût le mettrait

Toute reproduction sans autorisation du Centre français d’exploitation du droit de copie est strictement interdite.
© Techniques de l’Ingénieur, traité Génie nucléaire B 3 070 − 23
MÉTHODES DE CALCUL NEUTRONIQUE DE CŒUR _____________________________________________________________________________________________

hors de portée d’utilisation. Pour contourner cette difficulté, il faut

avoir recours à une méthodologie permettant d’apporter les correc- Annexe — Équation générale
tions nécessaires au calcul de base, approchant au mieux la réalité
de la configuration à traiter. À titre d’exemple, le tableau suivant
du transport
montre le remontage de la valeur de la réactivité du cœur de pre-
mière divergence de Superphénix et la figure 6 présente la modé-
lisation (x, y ) utilisée pour le traitement fin de l’hétérogénéité des Bilan neutronique
barres SCP (Système de Commande Principal) de Superphénix.
Afin d’établir l’équation du bilan neutronique, on cherche à
exprimer de façon générale l’évolution de la population des neu-
trons dans un milieu multiplicateur quelconque.
Les variables qui caractérisent l’état des neutrons à l’instant t
sont [3] [4] sa position r, sa vitesse v et sa direction Ω.
La population neutronique est donc représentée par une distri-
bution notée : n (r, v, Ω, t ) d3r dv d2Ω dt.
Il s’agit ici d’une densité angulaire de neutrons qui est une dis-
tribution par rapport aux variables r, v et Ω et une fonction du
temps [nous négligeons, par ailleurs, les neutrons retardés (§ 2)].
Pour exprimer la variation ∂n/∂t de la population neutronique n
en r entre t et t + dt, on établit le bilan des apparitions et des dis-
paritions de neutrons de vitesse v à dv près et de direction Ω à d2Ω
près dans l’élément de volume d3r entourant r :
Figure 6 – Traitement de l’hétérogénéité des barres
du SCP de Superphénix ∂
[n (r, v, Ω, t + dt ) – n (r, v, Ω, t )] d3r dv d2Ω = n(r, v, Ω, t)
∂t

(unité en 10–5) Termes du bilan

Résultat de calcul hexagonal 3 D en 25 groupes – 1 200
Les phénomènes physiques à prendre en compte pour établir ce
Correction de maillage : bilan dans l’élément de volume considéré sont explicités ci-des-
— radial ............................................................... + 147 sous.
— axial ................................................................. + 96
Correction de transport :
— radial ............................................................... + 366 Entrées-sorties
— axial ................................................................. + 423
Corrections d’hétérogénéité................................... – 52 Le bilan des entrées et des sorties de neutrons est représenté par
le terme appelé fuites.
Correction de la teneur en Pu 241
(au 07/09/1985) ........................................................ + 80 On somme ces fuites sur la surfaces dS entourant d3r :

« Formulaire »..........................................................
Valeurs toutes corrections effectuées ...................
+ 260
+ 120 (1)
Fuites =  dS
n ( r ,v , Ω , t )v Ω dS

(1) L’incertitude associée est de ± 310 × 10–5. qui, en utilisant le théorème d’Ostrogradsky [1], permet d’exprimer
la divergence du vecteur n (r, v, Ω, t ) v Ω par :

Pour certaines grandeurs, telles que puissances ponctuelles et div[n (r,v, Ω, t ) v Ω] = v Ω∇n (r, v, Ω, t )
taux de dommage, il n’y a pas en général de facteurs correctifs à où ∇n est le gradient de n (r, v, Ω, t ).
prendre en compte (dans les zones fissiles) : ceux-ci sont, au niveau
Si d3r est suffisamment petit pour que ∇n puisse y être
du projet, fondus dans les incertitudes. Elles sont donc issues du
considéré comme uniforme, alors le terme de fuites s’exprime par :
calcul de base.
Les incertitudes prises en compte dans le projet neutronique des Fuites = v Ω ∇n (r, v, Ω, t ) d3r dv d2Ω dt
RNR sont données à un intervalle de confiance de 2 σ et, en général,
se combinent quadratiquement. Les principales sont les incertitudes
sur les niveaux de réactivité, sur les taux de réaction (fission Pu239, Disparitions de neutrons
fission U235, fission U238, capture U238, etc.), sur les puissances
(puissances linéiques, puissances intégrées), sur l’effet Doppler et Toute collision d’un neutron de vitesse v et de direction Ω avec
les coefficients de réactivité, et sur l’effet de vide sodium. Une dis- un noyau de la matière contenue dans l’élément de volume d3r
cussion approfondie sur les incertitudes est présentée dans [81]. retournant r conduit :
La qualification du formulaire Carnaval-IV est basée sur l’inter- — soit à une absorption, le neutron disparaît alors totalement ;
prétation de dizaines d’expériences effectuées en France sur — soit à une modification de la direction ou de la vitesse du neu-
maquette critique et sur le retour d’expérience des essais de démar- tron (le plus souvent les deux à la fois), le neutron disparaît ainsi
rage et de redémarrage du réacteur Superphénix. On trouvera les de la population de vitesse v et de direction Ω.
principaux éléments de la validation et qualification de Carnaval-IV
dans [82] [83] [84].

Toute reproduction sans autorisation du Centre français d’exploitation du droit de copie est strictement interdite.
B 3 070 − 24 © Techniques de l’Ingénieur, traité Génie nucléaire
____________________________________________________________________________________________ MÉTHODES DE CALCUL NEUTRONIQUE DE CŒUR

Soit Σt(r, v, Ω, t ) la section efficace totale du milieu au point r et On préfère en général utiliser le flux neutronique, produit de la
à l’instant t, les disparitions de neutrons par collision correspondent densité angulaire de neutrons par leur vecteur vitesse, plutôt qu’une
à: distribution de population (car les taux de réactions s’expriment
Disparitions = Σt(r, v, Ω, t ) v n (r, v, Ω, t ) d3r dv d2Ω dt simplement par le produit d’un flux par une section efficace), soit :
Φ (r, v, Ω, t ) = vn (r, v, Ω, t )
On utilise plus couramment l’énergie E que la vitesse v du neu-
Arrivées par transfert en énergie tron. L’équation du bilan s’exprime alors par :
1 ∂Φ ( r, E, Ω, t ) 3 2
De tous les neutrons présents dans le volume d3r entourant r, de ---- ---------------------------------- = – Ω∇ Φ ( r, E, Ω, t )d rdEd Ω
v ∂t
vitesse v’ et de direction Ω’ quelconques, certains vont subir une
3 2
collision qui, par diffusion élastique ou inélastique, va les amener – Σ t ( r, E, Ω, t )Φ ( r, E, Ω, t )d rdEd Ω
à la vitesse v dans la direction Ω.

 
∞
La section efficace différentielle de transfert Σs(r, v’ → v, Ω’ → Ω, t ) 2
+ dE′ d Ω ′ Σ s ( r, E′ → E, Ω′ → Ω, t )x
caractérise ce phénomène. 0 4π
Par sommation sur toutes les vitesses et directions des neutrons
présents en r, on calcule ceux qui après choc vont avoir la vitesse 3 2
v et la direction Ω, soit : Φ′ ( r, E′, Ω′ , t ) d rdEd Ω
Transferts =

  + P ( r, E, Ω, t ) + Q ( r, E, Ω, t )
∞
2
dv′ d Ω′ Σ s (r, v′ → v, Ω′→Ω, t)
0 Cette équation, dite équation du transport ou équation de Boltz-
4π
mann, est différentielle par rapport au temps et à l’espace et inté-
3 2
grale par rapport à l’énergie et à la direction des neutrons.
v ′n′ ( r, v′, Ω′ , t ) d r dv d Ω dt
En régime stationnaire, productions et disparitions des neutrons
s’équilibrent, et le système reste stable. En l’absence de sources
extérieures, l’équation précédente peut alors s’écrire sous la forme
compacte :
Productions de neutrons A0 Φ 0 + P0 Φ 0 = 0
dans l’élément de volume
avec :
On distingue deux sortes de sources de neutrons.
A0 = –[Ω∇Φ0(r, E, Ω, t ) + Σ t (r, E, Ω, t ) Φ0 (r, E, Ω, t )] d3r dE d2Ω
Les sources de fission P (r, v, Ω, t ) dont l’intensité est proportion-

 
nelle au taux de réaction en chaque point et à chaque instant : ∞
2 3 2
Productions = P (r, v, Ω, t ) P0 = dE′ d Ω ′ Σ f ( r, E′ , t ) Φ0′ ( r, E′, Ω′ , t ) d r dE d Ω

 
0 4π
∞
2
= dv ′ d Ω ′ Σ f ( r, v ′ → v,  ′ →  , t )v′n′ ( r,v ′,  ′,t ) représentant respectivement les opérateurs intégraux de dispari-
0 4π tion (absorption plus fuites) et de production du système multipli-
cateur. On retrouve ainsi l’équation (10) dans le cas particulier
En général, ce terme se simplifie car on considère la production
λ0 = 1.
de neutrons de fission comme isotrope et on ne tient pas compte 1 ∂Φ ( r, E, Ω, t )
de la dépendance énergétique des neutrons incidents dans le spec- Si l’on définit :  = ---- ---------------------------------- – Q ( r, E, Ω, t )
v ∂t
tre de fission.
Les autres sources de neutrons, notées Q (r, v, Ω, t ), incluent, par en remplaçant dans l’équation du bilan les formes compactes pré-
exemple, les sources de démarrage. Dans un réacteur en fonction- cédentes, on retrouve aussi l’équation (1).
nement, ces sources sont en général faibles vis-à-vis de la produc- Cette équation est l’expression du bilan neutronique dans
tion de neutrons par les fissions ; elles sont le plus souvent l’approximation qui consiste à négliger les neutrons retardés. Dans
négligées lorsque l’on suppose le cœur critique. le cas plus général, la définition de  contient un terme relevant de
la présence de ces neutrons retardés. Comme il est montré en [17],
ce terme est extrêmement important pour les études cinétiques.
Bilan global

Le bilan des apparitions et des disparitions de neutrons au point

r s’écrit alors :
∂ 3 2
n ( r, v, Ω, t ) = – v Ω∇n ( r, v, Ω, t )d r dv d Ω
∂t
3 2
– Σ t ( r, v, Ω, t ) vn ( r, v, Ω, t )d r dv d Ω

 
∞
2
+ dv ′ d Ω ′ Σ s ( r, v′ → v, Ω′ → Ω, t )x
0 4π

3 2
v ′n′ ( r, v ′, Ω′ , t ) d r dv d Ω

+ P ( r, v, Ω, t ) + Q ( r, v, Ω, t )

Toute reproduction sans autorisation du Centre français d’exploitation du droit de copie est strictement interdite.
© Techniques de l’Ingénieur, traité Génie nucléaire B 3 070 − 25
 P
O
U
Méthodes de calcul neutronique R
de cœur
E
N
par Giovanni B. BRUNA
Docteur en Physique Nucléaire, Université de Gênes
Chargé de missions au Département Performances de Cœur de la société Framatome
S
et Bernard GUESDON
Ingénieur de l’École Nationale Supérieure de Mécanique et du Génie Atomique
A
Adjoint technique au Chef de Division Procédé de la société Framatome
V
Références bibliographiques O
Ouvrages généraux
[1] FRIEDMAN (B.). – Principles and techniques
Advances in nuclear science and technology,
vol. 9, Plenum Press, Londres et New York,
rements. In Proceedings of the Joint Interna-
tional Conference on Mathematical Methods
I
[2]
of applied mathematics. Dover Publications
Inc., New York, (1990).
ROZON (D.). – Introduction à la cinétique des
[15]
(1976).
GANDINI (A.). – Generalized perturbation
theory (GPT) methods. A heuristic approach. [23]
and Supercomputing in Nuclear Applica-
tions, Karlsruhe, (1993).
KATO (Y.), TAKEDA (T.) et TAKEDA (S.). –
R
réacteurs nucléaires. Éditions de l’École Poly- In J. Lewins, M. Becker Eds., Advances in Coarse-mesh correction of finite-difference
technique de Montréal, Montréal, (1992). nuclear science and technology, vol. 19, Ple- method for neutron diffusion calculations.
num Press, Londres et New York, (1987). Nucl. Sci. Eng. 61, p. 127-141, (1976).
[3] BUSSAC (J.) et REUSS (P.). – Traité de neutro-

[4]
nique. Hermann Ed, Paris, (1978).
RONEN (Y.). – CRC Handbook of nuclear reac-
[16] GOMIT (J.M.) et PLANCHARD (J.) Eds. – Série
des Bulletins de la Direction des Études et
Recherches EDF sur les méthodes de pertur-
[24] BRUNA (G.B.), CORNILUS (P.L.), POINOT (C.)
et VAN FRANK (C.). – Analysis of the main
physical & numerical aspects in the model-
P
[5]
tor calculations (3 volumes), CRC Press, Boca
Raton, Floride, (1986).
PLANCHARD (J.). – Méthodes mathéma-
bation en physique des réacteurs : série A,
no 1, (1985), série A, no 2, (1990), série C,
no 4, (1990), série A, no 1, (1991).
ling of control clusters for project calcula-
tions. In Proceedings of the IAEA Specialist
Meeting on Advanced Calculational Methods
L
tiques en neutronique. Collection des Études
et Recherches d’EDF, Eyrolles Ed., Paris,
(1995).
Démarche de conception
[17] BRUNA (G.B.) et SARGENI (A.). – A CPT clas-
for Power Reactors, Cadarache, sept. 1990.
Calculs du transport des neutrons
U
[6] WEINBERG (A.M.) et WIGNER (E.P.). – The
physical theory of neutron chain reac-
[Link] of Chicago Press, Chicago,
sical perturbation theory approach to slow
varying, small reactivity reactor kinetics. Pro-
gress in Nuclear Energy, vol. 28, no 2, p. 115-
[25] CARLSON (B.G.) et LATHROP (K.D.). – Trans-
port theory method of discrete ordinates. LA-
3251-MS (rev. Oct. 22, 1965) Los Alamos
S
(1968). 128, (1994). National Laboratory, Los Alamos, N.M,
(1965).
[7] ASH (M.). – Nuclear reactor kinetics (second Techniques de calcul
Edition). Mc Graw Hill Series in Nuclear Engi- [26] CARLSON (B.G.) et al. – Solution of the trans-
[18] MITA (J.R.) et BANASIAK (J.). – Derivation of
neering, Mc Graw Hill, (1979). port equation by the SN method. In Procee-
the neutron diffusion equation. In Y. Ronen et
dings of the U.N. International Conference on
[8] SPANIER (J.) et GELBARD (E.M.). – Monte E. Elias, Eds., Reactor physics and reactor
Peaceful Uses of Atomic Energy, Genève,
Carlo principales and transport problems. computations, Ben Gurion University of the
(1958).
Addison-Wesley Reading. Mass. (1969). Negev Press, (1994).
[27] KAVENOKY (A.). – Calcul et utilisation des
[9] LEWIS (E.E.) et MILLER Jr. (W.F.). – Computa- [19] SANCHEZ (R.), MONDOT (J.), STANKOVSKI
probabilités de première collision pour des
tional methods of neutron transport. A Wiley- (Z.), COSSIC (A.) et ZMIJAREVIC (I.). – A user
milieux hétérogènes à une dimension : les
Interscience Publication, J. Wiley & Sons, oriented, portable, modular code for multi-
programmes Alcoll et Cortina, CEA n-1077,
New York, (1978). group transport assembly calculations. In
2 - 1996

(1969).
[10] MARCHOUK (G.) et SHAYDOUROV (V.). – Raf- Proceedings of the International Topical Mee-
ting on Advances in Reactor Physics, Mathe- [28] BRUN (A.M.). – Généralisation de la méthode
finement des solutions des schémas aux dif-
matics and Computation, Paris, (1987). des probabilités de collisions pour tenir
férences. Éditions de Moscou (Trad.
compte soit du gradient du flux soit du choc
française), (1983). [20] COSTE (M.), SANCHEZ (R.), STANKOVSKI
linéairement anisotrope. Thèse, Université
[11] VARGA (R.S.). – Matrix iterative analysis. (Z.) et VAN DER GUCHT (C.). – Apollo II
d’Orsay, (1970).
Prentice-Hall Series in Automatic Computa- assembly spectrum code : new features. In
Proceedings of the International Topical Mee- [29] ENGLE Jr. (W.W.). – A user’s manual for
tion, Englewood Cliffs, N.J, (1962).
Doc. B 3 070

ting on Advances in Mathematics, Computa- ANISN. A one dimensional discrete ordinates

[12] WACHSPRESS (E.L.). – Iterative solution of transport code with anisotropic scattering.
tions and Reactor Physics, Pittsburgh, (1991).
elliptic systems and applications to the neu- USAEC Rep. K-1693, Union Carbide Corpora-
tron diffusion equations of reactors physics. [21] BRUNA (G.B.), CORNILUS (P.L.), GROSS-
tion, Oak-Ridge, Tenn, (1967).
Prentice-Hall Series in Automatic Computa- HANS (M.), NOBILE (M.) et VERGAIN (M.L.). –
Apollo II utilisation for project calculations. In [30] RHOADES (W.A.), SIMPSON (D.B.), CHILDS
tion, Englewood Cliffs, N.J, (1966).
Proceedings of the International Topical Mee- (R.L.) et ENGLE (W.W.). – The DOT IV, two
[13] NAKAMURA (S.). – Computational methods dimensional discrete ordinates transport
ting on Advances in Reactor Physics, Char-
in engineering and science. With applications code with space dependent mesh and qua-
lestons, 8-11 mars 1992.
to fluid dynamics and nuclear systems. A drature. ORNL/TM-6529, Oak-Ridge National
Wiley-Interscience Publication, J. Wiley & [22] RIMPAULT (G.), COLACIOPPO (V.) et
Laboratory, Oak-Ridge, Tenn., (1979).
Sons, New York, (1977). ROWLANDS (J.L.). – Recent methods in the
european cell code Ecco with applications to [31] RHOADES (W.A.) et CHILDS (R.L.). – Tort.
[14] GREENSPAN (E.). – Developments in pertur- Three dimensional Oak-Ridge discrete ordi-
the Sefor and Superphénix Doppler measu-
bation theory. In J. Lewins, M. Becker Eds.,

Toute reproduction sans autorisation du Centre français d’exploitation du droit de copie

est strictement interdite. − © Techniques de l’Ingénieur, traité Génie nucléaire Doc. B 3 070 − 1
P MÉTHODES DE CALCUL NEUTRONIQUE DE CŒUR _____________________________________________________________________________________________
O
U nates neutron/photon transport code. In Pro-
ceedings of the International Topical Meeting
on Advances in Mathematics, Computation
[48] WAGNER (M.R.) et al.. – A non linear exten-
sion of the nodal expansion method. In Pro-
ceedings of the International Topical Meeting
[65] PORTA (J.), DORIATH (J.Y.) et MARTINET (B.).
– Internal interface perturbations : calculatio-
nal method using the eigenvalue implicit

R [32]
and Reactor Physics, Pittsburgh, Penn.,
(1991).
HOFFMANN (A.), JANPIERRE (F.), KAVE-
on Advances in Mathematical Methods for
the Solution of Nuclear Engineering Pro-
blems, Munich, 27-29 avril 1981. [66]
expansion of the conditioned problem. Nucl.
Sci. Eng., 92, p. 367-371, (1986).
GANDINI (A.). – Implicit and explicit higher
NOKY (A.), LIVOLANT (M.) et LORAIN (H.). – [49] REMPE (K.R.), SMITH (K.S.) et HENRY (A.F.). – order perturbation methods for nuclear reac-
Apollo, code multigroupe de résolution de Simulate 3 pin power reconstruction : metho- tof analysis. Nucl. Sci. Eng., 67, p. 347-355,
l’équation du transport pour les neutrons dology and benchmarking. In Proceedings of (1978).

E [33]
thermiques et rapides. CEA N-1610, (1972).
SANCHEZ (R.) et al. – Apollo II. A user-orien-
the International Reactor Physics Confer-
ence, Jackson Hole, Wyoming, (1988).
[67] GANDINI (A.). – On the standard perturbation
theory. Nucl. Sci. Eng., 79, p. 426-439, (1981).
ted, portable, modular code for multigroup [50] HOBSON (G.H.) et AIGLE (R.). – Nodal code
N transport assembly calculation. In Procee-
dings of the IAEA Specialist Meeting on
developments at Gramatome/BWFC. In Pro-
ceedings of the International Topical Meeting
[68] GOMIT (M.), PLANCHARD (J.) et SARGENI
(A.). – La méthode des pseudo-
harmoniques : théorie et applications. Bulle-
Advanced Calculational Methods for Power on Advances in Reactor Physics, « Reactor tin de la Direction des Études et Recherches,
Reactors, Cadarache, sept. 1990. Physics Faces in the 21st Century », EDF, no 1, (1985).
[34] ASKEW (J.R.), FAYERS (F.J.) et KEMSHELL Knoxville, 11-15 avril 1994.
[69] PALMIOTTI (G.). – The FRENCH (Flux Recons-

S (P.B.). – A general description of the lattice

code WIMS. J. Br Nucl. Energy Soc., 5, 564,
(1966).
Méthodes de perturbation
[51] WIGNER (E.P.). – Effect of small perturbations
on pile period. Chicago Report CP6G-3048,
tructed Explicitly by Neutronic Calculated
Harmoniques) method in multidimensional
and multigroup theory. Ann. Nucl. Energy,

A [35]

[36]
DUBI (A.) et GERSTL (S.A.W.). – Nucl. Sci.
Eng, 76, 2, (1980).
TANG (J.S.), HOFFMAN (T.J.) et STEVENS
[52]
(1945).
USACHEV (L.N.). – Equation of the impor- [70]
no 4, p. 167-176, (1987).
BRUNA (G.B.) et SARGENI (A.). – A computa-
tional technique for evaluating eigenfunc-
tance of neutrons, reactor kinetics and the
V [37]
(P.N.). – Nucl. Sci. Eng., 64, 837, (1977).
NIMAL (J.C.). – Programme de Monte-Carlo
theory of perturbation. In Proceedings of the
First ICPUAE, vol. 5, United Nations, Genève,
tions of symmetrical nuclear systmes. Ann.
Nucl. Energy vol. 21, no 12, p. 745-758,
(1994).
polycinétique à trois dimensions Tripoli 01. (1955).
O [38]
Tomes 1 à 7, CEA N-1919, (1976).
MCMP. A general Monte-Carlo code for neu-
[53] PORTA (J.), DORIATH (J.Y.), MOURGUES (A.),
NOBILE (M.), VALLEE (A.) et BURNA (G.B.). –
[71] AIGLE (R.), BRUNA (G.B.) et SARGENI (A.). –
On the completeness of the multigroup
eigenfunctions set of a reactor system Boltz-
I tron and photon transport. LA-7396-M
(revised), Los Alamos National Laboratory,
Los Alamos, N.M., (1979).
Sensivity studies on high-conversion power
reactors, via perturbation methods and infor-
mation conditioning. Nucl. Sci. Eng., 95,
mann operator. Ann. Nucl. Energy, vol. 21,
no 8, p. 445-460, (1994).

R [39] BLOMQUIST (R.N.), LELL (R.M.) et GELBARD

(E.M.). – VIM. A continuous energy Monte-
Carlo code at ANL. ORL/RSIC-44, Oak-Ridge
[54]
p. 266-281, (1987).
USACHEV (L.N.). – Perturbation theory for
the breeding factor and other ratios of diffe-
[72] WILLERMOZ (G.), BRUNA (G.B.), CASTELLI
(R.) et BETHOUX (P.). – Studies on the sensi-
vity of PWR core parameters to fuel manufac-
National Laboratory, Oak-Ridge, Tenn., rent processes in a reactor. In Atomnaia Ener- turing uncertainties via statistical and
(1980). giya, 15, p. 472-181, (1963). perturbation methods. In Proceedings of the
ANS International Topical Meeting on Advan-
Méthodes de résolution [55] LEWINS (J.). – Importance, the adjoint func-

P [40] SCHMIDT (F.A.R.). – Finite element applica-

tion to global reactor analysis. In Procee- [56]
tion. Pergamon Press, (1965).
GANDINI (A.). – A generalized perturbation
[73]
ces in Reactor Physics, Knoxville, (1994).
MALDONADO (G.I.), TURINSKY (P.J.), KRO-
PACZEK (D.J.) et PARKS (G.T.). – Constrained
method for bi-linear functionals of the real
L dings of the International Topical Meeting on
Advances in Mathematical Methods for the
Solution of Nuclear Engineering Problems.
and adjoint neutron fluxes. Journal of
Nuclear Energy, 21, p. 755-765, (1967).
minimisation of fresh reload fuel enrichment
within the Formosa code’s nodal GPT and
optimization framework. In Proceedings of

U [41]
Munich, avril 1981.
DAVIERWALLA (D.M.), PELLONI (A.) et ALTI-
PARMOKOV (D.). – Finite elements in multidi-
[57] BRUNA (G.B.), AIGLE (R.) et LOUSSOUARN
(O.). – A generalized perturbation theory for
the reconstruction of the axial offset in a
the ANS International Topical Meeting on
Advances in Reactor Physics, Knoxville,
(1994).

S mensional nuclear diffusion. In Proceedings

of the International Topical Meeting...
Munich, (1981). [58]
pressurized water reactor. Nucl. Sci. Eng. 105,
p. 290-296, (1990).
GANDINI (A.). – The role of the control varia-
[74] MALDONADO (G.I.) et TURINSKY (P.J.). –
Application of nonlinear nodal diffusion
generalized perturbation theory to nuclear
[42] LAUTARD (J.J.). – New finite element repre- ble in the heuristically based generalized per- fuel reload optimization. Nucl. Tech., 110,
sentation for 3D reactor calculations. In Pro- turbation theory. In Proceedings of the ANS p. 198-219, (1995).
ceedings of the International Topical International Topical Meeting on Advances in
Meeting... Munich, (1981). Reactors Physics and Reactor Computation, Schémas de calcul
Portland, (1995). [75] VALLEE (A.), FRANCILLON (G.) et PELET (J.).
[43] BRUNA (G.B.) et SARGENI (A.). – A computa-
[59] MITANI (H.). – Higher order perturbation – SCIENCE : a new reactor physics codes pac-
tional technique for evaluating eigenfunc-
methods in reactor calculation. Nucl. Sci. kage approach. In Proceedings of the Topical
tions of symmetrical nuclear systems. Ann.
Eng. 51, p. 180-188, (1973). Meeting on Advances in Reactor Physics,
Nucl. Energy, vol. 21, no 12, p. 745-758,
Charleston, 8-11 mars 1992.
(1994). [60] KALLFELZ (J.M.), BRUNA (G.B.), PALMIOTTI
(G.) et SALVATORES (M.). – Burn-up calcula- [76] GIRIEUD (P.), LOUSSOUARN (O.), ROSHD
Méthodes nodales avancées tions with time-dependent generalized per- (M.) et VALLEE (A.). – An intelligent front-end
[44] KOEBKE (K.). – Advances in homogenization turbation theory. Nucl. Sci. Eng. 62, p. 304- for core computations. In Proceedings of the
and dehomogenization. In Proceedings of the 309, (1977). Topical Meeting on Advances in Reactor Phy-
International Topical Meeting on Advances in sics, Charleston, 8-11 mars 1992.
[61] CACUCI (D.G.), WEBER (C.F.), OBLOW (E.M.)
Mathematical Methods for the Solution of [77] GIRIEUD (P.). – SCIENCE : the new Frama-
et MARABLE (J.H.). – Sensivity theory for
Nuclear Engineering Problems, Munich, 27- tome 3D nuclear code package for safety ana-
general systems on nonlinear equations.
29 avril 1981. lysis. Presented at the ENC’94 Meeting, Lyon,
Nucl. Sci. Eng., 75, p. 88-110, (1980).
[45] SMITH (K.S.). – QPANDA : an advanced nodal oct. 1994.
[62] GANDINI (A.). – Generalized perturbation
method for LWR analysis. Trans. Am. Nucl. [78] DORIATH (J.Y.), McCALLIEN (C.W.), KIEFHA-
theory for nonlinear systems from the impor-
Soc., 50, 532, (1985). BER (E.), WEHMANN (U.) et REUNIER (J.M.).
tance conservation principle. Nucl. Sci. Eng.,
[46] LAWRENCE (R.D.). – Progress in nodal 77, p. 316-343, (1981). – ERANOS 1 : the advanced european system
methods for the solution of the neutron diffu- of codes for reactor physics calculations. In
[63] RAHNEMA (F.) et POMRANGING (G.C.). –
sion and transport equations. Progress in Proceedings of Joint International Confer-
Boundary perturbation theory for inhomoge-
Nuclear Energy, vol. 17, no 3, p. 271-301. ence on Mathematical Methods and Super-
neous transport equations. Nucl. Sci. Eng.,
(1986). computing in Nuclear Applications,
84, p. 313-319, (1983).
[47] SMITH (K.S.). – Assembly homogenization Karlsruhe, (1993).
[64] RANHEMA (F.). – Internal interface perturba-
techniques for light water reactor analysis. [79] RIMPAULT (G.), COLACIOPPO (V.) et
tions in neutron transport theory, Nucl. Sci.
Progress in Nuclear Energy, vol. 17, no 3, ROWLANDS (J.L.). – Recent methods in the
Eng. 86, p. 76-90, (1984).
p. 303-335, (1986). european code Ecco with applications to the
Sefor and Superphenix doppler measure-

Toute reproduction sans autorisation du Centre français d’exploitation du droit de copie

Doc. B 3 070 − 2 est strictement interdite. − © Techniques de l’Ingénieur, traité Génie nucléaire
____________________________________________________________________________________________ MÉTHODES DE CALCUL NEUTRONIQUE DE CŒUR
P
O
ments. In Proceedings of Joint International
Conference on Mathematical Methods and
Supercomputing in Nuclear Applications,
[83] BERGEONNEAU (Ph.) et al. – Analysis and
interpretation of the neutronic experimental
results at the Superphenix re-start-up in
[86] CHAUVIN (J.P.), GRANGET (G.), MARTINI
(M.), MONDOT (J.), LEFEBVRE (J.C.) et VAL-
LEE (A.). – Validation of the experimental
U
[80]
Karlsruhe, (1993).
GREENSPAN (E.). – A generalized source
multiplication method for determining reacti- [84]
1989. In Proceedings of the Physor 90 Mee-
ting, Marseille, (1990).
CABRILLAT (J.C.) et al. – Comparison of cal-
methods used in the Epicure program and
their associated uncertainties. In Proceedings
of the International Topical Meeting on
R
vity. Nucl. Sci. Eng., 55, p. 10-103, (1975). culations and measurments of reactions Advances in Reactor Physics, Charleston,
rates in the outer regions of Superphenix. In (1992).
[81] BERGEONNEAU (Ph.) et al. – R and D pro-
gramme in France in support of design stu- Proceedings of the Physor 90 Meeting, Mar- [87] BRUNA (G.B.). – Validation of pin power cal-
dies devoted to new FBR’s core concepts. In
Proceedings of the Physor 90 Meeting, Mar- [85]
seille, (1990).
MONDOT (J.), GAUTIER (J.C.), CHAUCHE-
culations on an experimental basis. In Pro-
ceedings of Workstop on Nuclear Reactors-
Physics, Design and Safety, Trieste, (1995).
E
seille, (1990). PRAT (P.), CHAUVIN (J.P.), GARZENNE (C.) et
[82] GAUTHIER (J.C.) et al. – The Balzac Program :
summary of the comparison between the
LEFEBVRE (J.L.). – EPICURE : an experimen-
tal program devoted to the validation of the
[88] BRUNA (G.B.), DOUCET (M.), NOBILE (M.) et
VALLEE (A.). – Up-dating and validation of a
N
main calculated and exprimental results. In calculational schemes for plutonium recy- project MOX computational chain. In Procee-
Proceedings of the Physor 90 Meeting, Mar- cling in PWRs. In Proceedings of the Physor dings of the 1992 Topical Meeting on Advan-
seille, (1990). 90 Meeting, Marseille, (1990). ces in Reactor Physics, Charleston, 8-11 mars
1992.
S
A
V
O
I
R

P
L
U
S

Toute reproduction sans autorisation du Centre français d’exploitation du droit de copie

est strictement interdite. − © Techniques de l’Ingénieur, traité Génie nucléaire Doc. B 3 070 − 3