Study of sensitivities of Monte Carlo simulations to nuclear resonance parameters

IN

DEGREE PROJECT

ENGINEERING PHYSICS,

SECOND CYCLE, 30 CREDITS

,

STOCKHOLM SWEDEN 2018

Study of sensitivities of Monte

Carlo simulations to nuclear

resonance parameters

ELIAS VANDERMEERSCH

Study of sensitivities of Monte

Carlo simulations to nuclear

resonance parameters

VANDERMEERSCH Elias

2018-05-17

Master’s

Thesis

Academic advisers

Pierre Tamagno (CEA)

Jan Dufek (KTH)

Adrien Bidaud (PHELMA)

Abstract | 1

Abstract

The nuclear data used for Monte Carlo simulation of nuclear reactors are retrieved from evaluated nuclear libraries. These data, and particularly the cross sections, are regularly corrected by integral experiments assimilation. Nevertheless, cross sections found in these libraries are computed from model parameters which describe interactions between the compound nucleus and neutrons. It makes then more sense to study integral experiment Monte Carlo simulation sensitivities to these model parameters. As cross sections evaluations depend on the energy domain studied, the work presented here is restrained to the Resolved Resonance Range, where the R-Matrix formalism can be applied.

The main purpose of this master’s thesis was to extend the TRIPOLI-4® _{Iterated Fission} Probability method ([1], [2]) developed by G. Truchet [3], to resonance parameters perturbations. This method uses the adjoint flux to compute Collision-based Exact Perturbations and sensitivities to cross sections, geometry and concentrations.

The first step was to compute sensitivities by direct finite difference, using critical Monte Carlo simulations, to obtain reference sensitivities. Such computations require to modify the libraries files used by TRIPOLI-4®_{. Two simulations are done, with one resonance parameters slightly modified.} Using these modified libraries, the exact perturbations between modified libraries are obtained, and the sensitivities to resonance parameters can be calculated. This method has the advantage of converging on the reactivity variation, but is laborious, and needs new Monte Carlo simulations for each studied parameter.

A coupling between CONRAD derivation functions and the IFP method post-processing code have then be implemented, to directly compute cross sections sensitivities to resonance parameters. This improvement offers the possibility to calculate with only one Monte Carlo simulation, sensitivities to all the resonance parameters of all the nucleus present in the system studied. Multiple tests have been conducted to ensure the reliability of the values obtained with 239_Pu, 240_Pu, 235_{U and} 238_{U. It appears that, with a proper number of direct batch used, and IFP simulation} parameters employed, results obtained are consistent with the direct perturbation method and converged faster, which proves the potential of this technique once correctly set up. This work also indicates possible upgrades and studies to conduct, based on the results obtained.

Keywords

Résumé | 3

Résumé

La simulation de réacteurs par code Monte Carlo nécessite l’utilisation de données nucléaires, généralement fournies sous forme de bibliothèque. Ces données, en particulier les sections efficaces, sont régulièrement corrigées par assimilation d’expérience intégrales. Néanmoins, ces sections efficaces sont elles-mêmes évaluées à partir de différents paramètres de modèle, qui caractérisent les interactions nucléaires susceptibles de se produire. Il est alors physiquement plus juste d’utiliser ces paramètres, et donc d’étudier la sensibilité d’un réacteur à ces derniers.

Selon le domaine d’énergie étudié, plusieurs formalismes permettent de construire les sections efficaces à partir de ces paramètres, dont celui de la matrice R, couramment utilisé pour reconstruire les sections efficaces dans domaine des résonances résolues. L’ensemble de ce travail se focalisera sur ce domaine.

L’objectif principal de ce stage était d’étendre la méthode des IFP ([1], [2]) développée par G. Truchet [3] au calcul de sensibilités aux paramètres de résonances.

I. Calculs de sensibilités

Trois méthodes furent utilisées durant ce stage pour calculer des sensibilités aux paramètres de résonance :

1. La méthode la plus évidente est de réaliser deux simulations Monte Carlo, l’une dans un état de référence, et l’autre dans un état perturbé, puis de les différencier. Pour un paramètre de résonance 𝛤𝑖, et deux états ref et pert :

𝜕𝜌 𝜕Γ𝑖

≈ 𝜌𝑟𝑒𝑓− 𝜌𝑝𝑒𝑟𝑡 Γ𝑟𝑒𝑓− Γ𝑝𝑒𝑟𝑡

(1) Cette technique permet d’obtenir un résultat qui peut être considéré comme de référence pour de petites variations. Toutefois, elle a le désavantage de nécessiter une modification des bibliothèques utilisées, ainsi que la réalisation d’au mieux une simulation Monte Carlo pour chaque sensibilité à étudier. Elle fut utilisée comme référence dans ce travail, afin de vérifier le bon fonctionnement des différents codes exécutés et développés.

2. La seconde détourne la méthode IFP développée par G. Truchet: Cette méthode utilise le flux adjoint pour calculer les sensibilités de simulation Monte Carlo d’expériences intégrales à des perturbations de concentration, volume, section efficace … En effet, l’équation 2 nous montre que pour de petites perturbations, il suffit de connaitre le flux 𝜙𝑟 de l’état de référence, et le flux adjoint 𝜙𝑝† _{de l’état perturbé pour déterminer la} variation de réactivité engendrée par la perturbation. Il n’y a donc plus besoin de réaliser une simulation complète mais de multiples petites simulations, de stocker les données de collisions et de calculer les flux.

𝛿𝜌 = 〈𝜙𝑝

†_{, [𝜆𝛿𝑃 − 𝛿𝐾]𝜙} 𝑟]〉

Résumé | 4

Si cette méthode est appliquée sur des sections efficaces reconstruites à partir de paramètres de résonances modifiés, on obtient la perturbation exacte engendrée par la modification des paramètres, et donc la sensibilité recherchée. Contrairement à la méthode 1, où la convergence est faite sur la valeur de la réactivité de chaque état, c’est sur la variation de réactivité que le calcul converge.

Cette méthode nécessite, comme la précédente, l’utilisation de bibliothèques modifiées pour chaque paramètre étudié. Elle est ici utilisée pour vérifier si la méthode IFP permet d’obtenir des perturbations équivalentes à celles des calculs directs. Des erreurs persistantes ont permis de déceler un problème dans le code du calcul du flux adjoint du post-traitement JIMMY2, qui est utilisé pour calculer les sensibilités à partir de fichiers de collision générés par TRIPOLI-4®_{. Ce problème provenait d’une erreur} simple, issue de modifications du code ultérieures aux travaux de G. Truchet. Après correction, les résultats obtenus sont bien similaires au calcul direct.

3. La dernière méthode fut celle développée durant ce travail : Le but était d’implémenter le calcul des sensibilités des sections efficaces aux paramètres de résonance directement dans la méthode IFP, et de propager ces sensibilités dans celles calculées (sensibilités de la réactivité au sections efficaces) par la méthode 2. Pour un paramètre de résonance 𝛤, un isotope i, et chaque réaction possible r, la sensibilité au paramètre 𝛤𝑖 s’exprime par :

𝑆𝜌,Γ_i = ∫ Γi⋅ 𝜕𝜌(𝐸, Ω) 𝜕Γi 𝑑𝐸𝑑Ω 𝐸,Ω = ∫ ∑ 𝜕𝜌(𝐸, Ω) 𝜕𝜎𝑖,𝑟(𝐸, Ω)⋅ 𝜕𝜎𝑖,𝑟(𝐸, Ω) 𝜕Γi ⋅ Γi 𝑟 𝑑𝐸𝑑Ω 𝐸,Ω (3)

La méthode IFP calculant le terme 𝜕𝜌(𝐸,Ω)

𝜕𝜎𝑖,𝑟(𝐸,Ω)⋅ 𝜎𝑖,𝑟(𝐸, Ω), on peut remonter à la sensibilité recherchée en implémentant le terme 𝜕𝜎𝑖,𝑟(𝐸,Ω)

𝜕Γ_i ⋅ Γi

𝜎_𝑖,𝑟(𝐸,Ω). Ce calcul est réalisé par des fonctions de CONRAD, couplées au code de post-traitement JIMMY2.

Pour calculer des sensibilités de référence, il a fallu générer des fichiers d’évaluation modifiés. Le code CONRAD possède un parseur de fichiers ENDF qui reconstruit les trois premières MF à partir d’un fichier d’évaluation de référence, et d’un fichier théorique (qui recense l’intégralité des paramètres nucléaires correspondant à un formalisme défini, tel Reich-Moore). Ce parseur a été modifié pour changer la MF2 d’un fichier ENDF, en prenant soin de corriger les différents liens inter-MF, de supprimer les fichiers concernant les covariances (non-utilisés ici, et qui ralentissent le traitement), et d’en d’indiquer par des flags la nécessité de reconstruire les sections efficaces à partir des paramètres de la MF2. La cohérence du fichier reconstruit est vérifiée par les outils Checkr et Stanef de l’IAEA.

Résumé | 5

II. Implémentation du calcul des sensibilités aux paramètres de résonance

Différentes classes de CONRAD permettent le calcul des sensibilités des sections efficaces aux paramètres de résonance. Ces classes furent liées au post-traitement JIMMY2, utilisé dans la méthode IFP. Le calcul des sensibilités est effectué de la façon suivante :

• Reconstruction des sections étudiées à 0K, en générant une grille d’énergies qui permet d’obtenir les sections interpolables sur cette grille à une précision prédéfinie.

• Calcul des dérivées à 0K.

• Elargissement des sections par prise en compte de l’effet Doppler, ajustement de la grille d’énergies.

• Calcul des dérivées à 300K.

• Pondération des sensibilités aux sections par les dérivées aux paramètres à 300K (JIMMY2 calcule les sensibilités à chaque collision (à une énergie E), pour chaque type de réaction).

Quelques optimisations ont été mises en place afin d’accélérer le calcul total. En effet, la génération de la grille consomme énormément en temps et en ressource : Pour chaque point, il faut reconstruire la section, calculer les dérivées de chaque paramètre, vérifier si le pas est suffisant, etc... Dans le cas de l’étude de la sensibilité à un seul paramètre, l’utilisateur peut indiquer en entrée le paramètre désiré, ainsi que le domaine de résonances sur lequel l’influence de ce paramètre est étudiée. Cette restriction permet de réduire grandement la quantité de calcul effectuée sans altérer le résultat. En effet, l’impact de la perturbation d’un paramètre sera nul sur une résonance éloignée. Pour la première résonnance du 239_{Pu, la variation de la sensibilité calculée (entre un calcul sur une} résonance voisine et l’ensemble du RRR) est inférieure au centième de pcm. En revanche, la quantité de mémoire utilisée est fortement réduite (de moins d’un GB pour un domaine de 10 résonances à plus d’une centaine pour 1000). L’utilisation de fichiers stocks allégés permet de même de réduire considérablement la quantité de mémoire nécessaire (contrairement à plus de 240GB pour des fichiers stocks). L’utilisateur a toujours la possibilité de réaliser un calcul de sensibilité/perturbation aux sections efficaces dans la version modifiée durant ce stage.

À noter que ce développement n’est pas optimum : En effet, les dérivées à tous les paramètres de résonances du domaine sont calculés à chaque exécution. La réduction du domaine permet donc de réduire la taille de la grille en énergie, et surtout le nombre de dérivées calculées. L’impact sur les ressources consommées par ces calculs inutiles étant relativement négligeable pour les petits domaines de résonance utilisés, l’implémentation n’a pas été plus optimisée durant ce stage. De même, un script bash est nécessaire pour calculer les sensibilités à plusieurs paramètres : les sensibilités aux sections sont donc recalculées pour chaque paramètre étudié. Une amélioration possible consisterait donc à ne calculer les dérivées qu’aux paramètres demandés, et de permettre d’indiquer plusieurs paramètres à étudier.

Résumé | 6

III. Vérification du fonctionnement et tests

Il est nécessaire d’identifier les paramètres optimaux à utiliser pour la méthode IFP. De nombreuses simulations ont été effectuées pour déterminer le nombre de batches direct, de générations par cycle adjoint, de neutrons, d’itérations et la valeur de la roulette russe à utiliser pour obtenir des résultats proches des calculs directs. Il faut néanmoins noter qu’il existe des corrélations entre les paramètres : la diminution du seuil roulette russe aura un effet similaire à augmenter le nombre de neutron par batches. Il n’existe donc pas une seule combinaison de paramètres optimaux. Les résultats obtenus après cette optimisation permettent de trouver des résultats très proches de ceux de référence. Les résultats comparés obtenus avec l’235_{U, pour une solution d’}235_{U et d’}238_{U sont} détaillés dans la Table 1. Il en ressort qu’en plus d’être relativement fiable, la méthode IFP permet l’obtention de valeurs bien plus précises que les calculs directs, pour un temps d’exécution inférieur.

Table 1 - Sensibilités à différents paramètres de résonnance d'une solution d'uranium

Différence finie (±1σ. pcm/%) IFP (±1σ. pcm/%) 𝐸𝜆2 37 ± 3 39.1 ± 0.3 𝐸𝜆3 -26 ± 6 -27.2 ± 0.5 𝐸𝜆7 -6 ± 6 -6.3 ± 0.6 Γ2 -15 ± 6 -14 ± 0.3 Γ3 17 ± 6 11.0 ± 0.1 Γ7 6 ± 6 2.91 ± 0.03

Ce travail a donc permis de démontrer la possibilité de coupler CONRAD à des méthodes Monte Carlo. On peut donc imaginer pouvoir, à l’avenir, simuler directement un réacteur à partir des paramètres décrivant les interactions nucléaires pour chaque domaine d’énergie (RRR, URR et Continuum), au lieu des sections efficaces. Un des avantages serait de pouvoir réduire les incertitudes liées à l’utilisation des fichiers ENDF. Il reste à contrôler la méthode précédente sur plus d’isotopes, en particulier des atomes légers. Le formalisme multilevel, bien qu’implémenté, n’a pas pu être testé. Il demeure donc ce travail à effectuer, même si ce formalisme tend à être abandonné au profit de la Matrice R réduite et de celui dit de Reich-Moore.

Mots-clés

Acknowledgments | 7

Acknowledgments

I would first like to thank my master’s thesis originator and supervisor, Pierre Tamagno, for the freedom he gave me, his help to solve unexpected bugs, and his assistance for this report. I also wish to express my gratitude to Jan Dufek and Adrien Bidaud, who have accepted to oversee this work from distance.

I am grateful to Cyrille De Saint Jean for allowing me to work at the LEPh, and to all the LEPh team for their goodwill, particularly Pascal Archier and Pierre Leconte for the tools and data they brought me.

A thousand thanks to all the students I met during these six months, starting with Axel: it was a pleasure to share this office with you, and I wish you all the best for your PhD. Good luck to Bastien, Paul, and Virginie, and thanks to Martin, Sylvain, Jordan, Daniela, Giorgio, Lino, Maxime, Timothée, for their friendship and time shared around a coffee. My regards to Ludovic, who kindly accepted to share some CPU during rush periods.

Table of contents | 9

Table of contents

Abstract ... 1

Keywords ... 1

Résumé ... 3

Mots-clés ... 6

Acknowledgments... 7

Table of contents ... 9

List of Figures ... 11

List of Tables ... 13

List of acronyms and abbreviations ... 15

Introduction ... 17

1

Evaluation of cross sections ... 19

1.1

Cross sections ... 19

1.1.1

Compound nucleus model ... 19

1.1.2

Resolved Resonances Range ... 21

1.1.3

Unresolved Resonances Range & Continuum ... 22

1.2

Evaluation of Nuclear Data ... 23

1.2.1

Evaluated Nuclear Data File ... 23

1.2.2

CONRAD ... 25

1.3

Summary ... 25

2

Exact perturbations in Monte Carlo simulation ... 27

2.1

Monte Carlo simulations ... 27

2.1.1

Principle ... 27

2.1.2

TRIPOLI-4

_{... 27}

2.2

Adjoint Flux ... 28

2.2.1

Definition ... 28

2.2.2

Iterated Fission Probability ... 29

2.3

Reactivity sensitivities ... 30

2.3.1

Eigenvalue difference method ... 30

2.3.2

Exact and first order Perturbation Theory ... 31

2.3.3

Collision based Exact Perturbation (CEP) ... 31

2.3.4

Sensitivities in JIMMY2 ... 34

2.4

Summary ... 35

3

Methodology ... 37

3.1

Sensitivities to resonance parameters ... 37

3.2

Approach ... 37

3.2.1

XDR files creation ... 38

3.2.2

Reference Calculation ... 38

3.2.3

Original IFP method testing ... 39

3.2.4

New IFP method ... 40

3.3

Environment ... 41

3.3.1

Libraries ... 41

Table of contents | 10

3.3.3

Hardware ... 42

3.4

Schedule ... 42

3.5

Summary ... 43

4

Analysis... 45

4.1

Direct calculation ... 45

4.2

Original IFP method results ... 47

4.3

Implemented method validation ... 47

4.3.1

Optimal simulation parameters ... 47

4.3.2

Performance and efficiency ... 51

4.3.3

_{Pu parameters PST sensitivities ... 51}

4.3.4

_{Pu parameters PST sensitivities ... 53}

4.3.5

_{U parameters “UST” (Uranium -Solution-Thermal)}

sensitivities ... 53

4.3.6

_{U parameters “UST” sensitivities ... 55}

Conclusion... 59

References ... 61

Appendix A: Results obtained for 𝚪𝑵 ... 63

Appendix B: Code developments ... 67

ENDF Files reconstruction ... 67

CONRAD “coupling” in the JIMMY2 post-processing code ... 68

List of Figures | 11

List of Figures

Figure 1.1 -

_{Pu total cross section from the JEFF3.2 evaluation. ... 19}

Figure 1.2- Energy states of the compound nucleus model. Elastic and

Inelastic scattering are depicted on the left [7]. ... 20

Figure 1.3 - Schematic of the nuclear data path from experimental and/or

theoretical sources to industrial applications, [9]. ... 23

Figure 2.1-Calculus of the reactivity change induced by all the possible

scattering reactions at the collision i. ... 33

Figure 2.2 - Calculus of the reactivity sensitivity to the

scattering reactions at the collision i. ... 34

Figure 3.1 - Process of XDR files generation ... 38

Figure 3.2 - Sensitivity obtained by direct Monte Carlo simulation ... 39

Figure 3.3 - Sensitivity obtained by the original IFP method ... 39

Figure 3.4 - New IFP Method dataflow ... 40

Figure 3.5 – Geometry of the PST used ... 41

Figure 3.6 - Initial planning of the project ... 42

Figure 3.7 - Real progress of the project. ... 43

Figure 4.1 - 𝐸𝜆1 variation impact on

_{Pu fission and capture cross ... 45}

Figure 4.2 - 𝛤𝑛1 variation impact on

_{Pu fission and capture cross ...46}

Figure 4.3 - Plots for various perturbed resonance parameters of the

reaction rate ratio ...46

Figure 4.4 - Impact of the direct batches number. ...49

Figure 4.5-

_{Pu cross sections sensitivities, and 𝐸𝜆 sensitivities for}

resonances 2 to 10, for 2000 direct batches. ... 52

Figure 4.6 -

_{Pu cross sections sensitivities, and sensitivities for each 𝐸𝜆,}

for 4000 direct batches. ... 52

Figure 4.7 -

_{Pu Cross sections sensitivities, and sensitivities for each 𝐸𝜆,}

for 2000 direct batches. ... 53

Figure 4.8 -

_{U 𝐸𝜆 sensitivities, for 8000 direct batches. ... 54}

Figure 4.9 - Zoom on the 7

_{and 10}

_{U resonance. ... 55}

Figure 4.10

_{U 𝐸𝜆 sensitivities, for 8000 direct batches. ... 56}

Figure 4.11 -. Zoom on some

_{U resonances, for 8000 direct batches. ... 56}

List of Tables | 13

List of Tables

Table 1.1 - Extract of various MF [10] ... 24

Table 1.2- Extract of various MT [10] ... 24

Table 3.1 - Calculators characteristics ... 42

Table 4.1 – PST-01-01 sensitivities to

_{Pu first resonance parameters 𝐸𝜆}

and 𝛤𝑁, computed by direct simulations ... 45

Table 4.2 - Some results obtained during testing of the IFP method.. ... 47

Table 4.3 - Relative difference (%) between direct and IFP computed

sensitivity to 𝐸𝜆1, for various batch size and number of

simulation. Time is normalized to the first column. ... 48

Table 4.4 - Relative difference (%) between direct and IFP computed

sensitivity to 𝐸𝜆1, for various Russian roulette (RR) threshold

and number of simulation. ...49

Table 4.5 - Relative difference (%) between direct and IFP computed

sensitivity to 𝐸𝜆1, for various IFP cycle length and number of

simulation. ... 50

Table 4.6 – PST-01-01 sensitivities to

_{Pu parameters. ... 50}

Table 4.7 - Effects of the resonance domain chosen and the type of file

used on the 𝐸𝜆1 sensitivity. ... 51

Table 4.8 – PST-01-01 sensitivities to

_{Pu parameters, for 2000 IFP}

direct batches. ... 53

Table 4.9 - UST sensitivities to

_{U parameters. ... 54}

Table 4.10 - UST sensitivities to

_{U parameters, for 8000 IFP direct}

List of acronyms and abbreviations | 15

List of acronyms and abbreviations

BBEP Black Body Exact Perturbations

CEA French Alternative Energies and Atomic Energy Commission CEP Collisions-based Exact Perturbations

CONRAD COde for Nuclear Reaction Analysis and Data assimilation

EDF Electricité de France

ENDF Evaluated Nuclear Data File GELINA Geel electron linear accelerator IAEA International Atomic Energy Agency

ICSBEP International Criticality Safety Benchmark Evaluation Project IFP Iterated Fission Probability

JEFF Joint Evaluated Fission and Fusion File JENDL Japanese Evaluated Nuclear Data Library LEPh Laboratory for Physical Studies

MF ENDF Evaluation Files

MT Reaction Nomenclature

PST Plutonium-Solution-Thermal

RRR Resolved Resonances Range

SERMA Service for Reactor Studies and Applied Mathematics

URR Unresolved Resonances Range

UST Uranium-Solution-Thermal

Introduction | 17

Introduction

The current issues of climate change and energetic needs [4] feed the debate about nuclear energy : it seems to be a worthy way to produce decarbonized electricity, but suffers from a public fear, mostly correlated to the consequences of previous nuclear accidents. Nuclear safety has been largely improved with the use of simulation codes, for conception and operation of a reactor. The importance of state-of-the-art code is regularly quoted by nuclear actors, as in GEN-IV forum [5] or by the IAEA.

According to the needs of nuclear safety, the French Alternative Energies and Atomic Energy Commission (CEA) continuously works to improve its nuclear codes. For these purposes, the Laboratory for Physicals Studies (LEPh) helps to integrate the feedback of integral experiments in nuclear data. Two kinds of approaches are currently used in reactor physics to simulate the neutron population:

Firstly, the deterministic approach (APOLLO-3), which uses approximations to solve the transport equation. Secondly, the stochastic approach (TRIPOLI-4®_{), which tracks individual} neutron histories by applying fundamentals nuclear physic laws. In both cases, to characterize the interaction of neutrons with matter, solvers use so-called nuclear data (angular distribution, spectra, emission, and particularly microscopic cross sections). These cross sections can be found in various databases (called evaluations), developed by multiple partnerships around the world. The values of these cross sections are not exact, as these are built from experimental data and theoretical models. Nevertheless, they are the image of the best knowledge we currently have.

Two techniques are used to evaluate nuclear data: Direct measurement, which requires specific experiments, as the Geel Electron Linear Accelerator (GELINA) [6], or assimilation of integral experiments. Adjustments on nuclear data are made, for example, by using data collected from nuclear reactor (like the reactivity).

Currently, assimilation of integral experiments is applied on cross sections energy groups. Such a method is practical but possesses multiple drawbacks: The multigroup structure is fixed, and then potentially limits the domain of application in which the evaluated data can be used. On the other hand, cross sections are partly evaluated from nuclear model parameters. Using these parameters in the assimilation process can offer more flexibility: Continuous cross sections can be obtained from adjusted nuclear parameters, which remove the limitation of the multigroup energy structure. It would make more sense to work with these parameters, and even overpass the evaluations file for any calculus on the long term.

A first step for the assimilation of integral experiences on resonance parameters is to calculate the sensitivities of integral parameters to resonance parameters. Such sensitivities can be computed with TRIPOLI-4®_{, but only by fastidious processes. It may be interesting to use current sensitivities} calculation methods implemented in TRIPOLI-4® _{to realise such computations, and thus offer an} efficient and quicker way to calculate the sensitivities of integral parameters to resonance parameters.

18 | Introduction

method can be used to evaluate nuclear data uncertainties propagation on integral experiments. Such a method generally need strong hypotheses. On the other hand, Monte Carlo sampling (stochastic approach) gives the possibility to overpass the deterministic hypotheses and to obtain reference sensitivities of integral values to nuclear data.

The purpose of this project is to extend the work of G. Truchet [3] and calculate exact sensitivities to resonance parameters (E𝜆, Γ𝑁, ΓF ...) using a stochastic approach. This work is restricted to the calculus of the nucleus’s sensitivities on its resolved resonances domain, using the Reich-Moore approximation (details in Chapter 2). It is a first step in the study of resonance parameters sensitivities and can be extended to the unresolved resonances and continuum domains.

Different steps must be achieved to manage the project:

• Reconstruct cross sections evaluation files with modified resonance parameters • Evaluate the sensitivities by direct calculations and test the IFP method

• Implement the calculus of cross sections sensitivities to resonance parameters • Verify this implementation running.

This thesis is split into five chapters: the first two parts present the theories and previous works linked to this subject, and the last three focused on what was accomplished during this master’s thesis.

In the first chapter, the relation between resonance parameters and cross sections - and how these data are made available for the scientific community – are introduced. The second chapter highlights methods used to compute sensitivities in Monte Carlo simulations, with an emphasis on the perturbations and the Iterated Fission Probability (IFP) methods using adjoint-flux and collisions. These two parts are the basis on which resonance parameters sensitivities computations will be achieve. The different ways to calculate these sensitivities and how they were processed and tested is detailed in chapter three. The chapter four summarizes the code developments realized, for data pre-processing and IFP method extension. Finally, results obtained with previous and new methods are presented in chapter five, with some details concerning the new method optimisation.

Evaluation of cross sections | 19

1 Evaluation of cross sections

1.1 Cross sections

Multiple interactions can occur on a nucleus inside a nuclear reactor, after a collision with a particle: capture, fission, scattering… The focus will be made here on reactions concerning incident neutrons. In order to describe the probability of a reaction to happen, physicists use so-called cross sections. Two kinds of cross section are generally used: microscopic cross section 𝜎, which characterizes interactions for a single nucleus, and macroscopic cross section Σ = 𝑁 ⋅ 𝜎, where the N is nucleus density, characterizing a medium.

Total cross section of 239_{Pu, which describes the total probability of a neutron to interact with a} 239_{Pu nucleus, is depicted on Figure 1.1. As it can be seen, cross sections are highly dependent on the} incident neutron energy. Three domains can be identified: The Resolved Resonance Region (RRR), the Unresolved Resonance Region (URR) and the Continuum. In the first one, multiple resonances can be observed. The density of resonances increases rapidly with the incident energy, so that at one point, it becomes impossible to distinguish two different resonances. A distinction is made between the URR and the continuum as the cross section could be determined by different methods in each range. Besides, URR contains parameters ensuring consistency of modelling with the RRR. The following parts introduce some basic aspects of the resonances theory.

1.1.1 Compound nucleus model

The compound nucleus model was introduced by Bohr in 1936 to explain interactions with an incident neutron [7]. In this model, the interaction process is split in two steps: the nucleus will first reach an excited state, and then decay by emitting particles. The excited nucleus is then said to have an entrance and an exit channel.

𝑋𝑍𝐴 + 𝑛 →𝐴+1 ∗𝑍𝑋 → 𝑌 + 𝑖 (1.1)

Figure 1.1 - 239Pu total cross section from the JEFF3.2 evaluation.

20 | Evaluation of cross sections

A channel is defined by:

• the particles pair 𝛼, including their mass (A, Z, E*), charges, spins I and parity. • the relative orbital angular momentum 𝑙

• its channel spin 𝑠, resulting of the coupling between spins of the pair 𝛼 • its total angular momentum and parity 𝐽π

It should be noticed that an interaction can only happen if the entrance channel total angular momentum equals the exit one (in respect of the momentum and parity conservation). A minimal separation energy Sn is necessary for a particle to be ejected from the nucleus, as depicted in Figure 1.2. Over this limit, the energy spectrum of the nucleus state is quasi-stationary. Such continuum should normally not create resonance. Yet, for energy superior to Sn, some states of the compound nucleus possess properties similar to the stationary ones [8]. These “decaying states” have, as every state, a life time 𝜏𝑠. According to the Heisenberg’s uncertainty principle, this life time is equivalent to an energetical “width”, Γ𝑠 _(1.5).

Γ𝑠_{⋅ 𝜏}

𝑠≅ ℏ (1.2)

Figure 1.2- Energy states of the compound nucleus model. Elastic and Inelastic scattering are depicted on the left [7].

Incident Neutron Energy

Evaluation of cross sections | 21

This width could be seen as the decay probability of the excited compound nucleus to go through all possible channels [8]. It means that Γ𝑠 _{could be split into partial widths, each of them} corresponding to the decay probability of each possible channel. If, for a decay state at energy E𝜆, the allowed interactions are only scattering reactions or radiative captures, the total width Γλ could be defined as:

Γ𝜆= Γ𝜆𝛾+ Γ𝜆𝑛+ Γ𝜆𝑛′ (1.3)

With Γ𝜆𝛾 the total radiative width, Γ𝜆𝑛the neutron width and Γ𝜆𝑛′ the inelastic neutron width. An intuitive relation of dependency seems to link reaction cross section to partial resonance width. Taking the following radiative capture (1.4) as an example, it is reasonable to think that the capture cross section σ𝑐 (E𝜆) will be reliant on Γ𝜆𝛾: the probability of a radiative capture at an energy E𝜆 to happen will be related to the probability of the compound nucleus at an energy E𝜆 to decay by emitting a photon 𝛾.

𝑛 01 +23994𝑃𝑢→24094𝑃𝑢∗ →24094𝑃𝑢+ 𝛾 (1.4) Such dependencies could be proven by advanced mathematics, using the collision matrix and R-Matrix formalisms [9]. The transition between the entrance and the exit channel could be described by the so-called collision matrix U. Considering respectively 𝑥𝑐′ and 𝑦𝑐 the “amplitude” of the waves in the outgoing and incoming channel:

𝑥𝑐′= − ∑ 𝑈𝑐′𝑐𝑦𝑐 𝑐

(1.5) The cross section of the transition could then be calculated using the U matrix coefficient,

𝜎𝑐𝑐′ = 𝜋 𝑘𝑐2

⋅ 𝑔𝐽⋅ |𝑒2𝑖𝑤𝑐𝛿𝑐𝑐′− 𝑈𝑐𝑐′|² (1.6) With 𝑔𝐽 the spin factor, depending on the spins I1 and I2 of the nucleus and the incident particle, k the wave number of the entrance channel, and 𝑤𝑐 related to the coulomb phase shift. The microscopic cross sections are then calculated by summing the partial cross section (1.6) over the entrance and exit channels concerned.

𝜎𝑡𝑜𝑡′= ∑ ∑ ∑ 𝜋 𝑘𝑐2⋅ 𝑔𝐽⋅ |𝑒 2𝑖𝑤𝑐_𝛿 𝑐𝑐′− 𝑈𝑐𝑐′|² ∀𝑐′ 𝑐′ ∀ 𝑐 𝑐 𝐽 (1.7)

The next step is to express the collision matrix in terms of the resonance parameters. This development depends on the energy domain studied.

1.1.2 Resolved Resonances Range

22 | Evaluation of cross sections

• Non-relativistic quantum mechanics is sufficient to describe the compound nucleus. • Only two-particles channels are treated.

• Beyond a distance ac, the potential field (considered Coulombian) in a pair c only depends on the nucleus/particle distance. This is translated in Equation (1.8) by the boundary condition B and the penetrability factor 𝑃𝑐, being perfectly described by 𝑘𝑐 and ac for each orbital angular momentum 𝑙.

The U Matrix could be expressed in terms of three other matrixes: The L-Matrix and the Ω-Matrix, which characterize external region (beyond ac), and the R-Matrix, associated with internal region (within ac).

𝑈𝑐𝑐′= Ω𝑐√𝑃𝑐{[1 − 𝑅(𝐿 − 𝐵)]−1[1 − 𝑅(𝐿∗− 𝐵)]} √𝑃𝑐′Ω𝑐′ (1.8) The potential inside the internal region is too complex to solve the Schrödinger equation. Nevertheless, with some boundary conditions B, the internal wave function could be express as a linear combination of the compound nucleus eigenstates. The R-Matrix coefficients could then be defined by [10]: 𝑅𝑐𝑐′(𝐸) = ∑ 𝛾𝜆𝑐𝛾𝜆𝑐′ 𝐸𝜆− 𝐸 𝜆 _(1.9)

With 𝛾𝜆𝑐 the reduced width, linked to the partial width by [8]:

Γ𝜆𝑐= 2 𝑃𝑐 𝛾𝜆𝑐² (1.10)

The use of the R-Matrix formalism highlights the presence of interferences between states and channels. Different approximations are used to diminish the complexity of the U-matrix coefficients:

• Single-level Breit-Wigner [11]: All the interferences terms are neglected, so that only one level per channel is considered. This formalism was mostly used for isolated resonances but is now considered obsolete.

• Multi-level Breit-Wigner [12]: Interferences between channels are neglected. Such formalism is used for light nucleus but could lead to non-physical cross sections.

• Reich-Moore [13]: In this formalism, only the interferences between 𝛾 channels are neglected. This leads to a reform of the R-Matrix, without the 𝛾 channels. The NEA now advises to use this formalism by default for new evaluations [14].

This part being only an introduction to the resonance theory, detailed cross sections expressions and derivatives to resonance parameters will not be developed here. Nevertheless, such expressions have been retrieved and calculated during this project, using CONRAD (see Section 0).

1.1.3 Unresolved Resonances Range & Continuum

Evaluation of cross sections | 23

transmission coefficients [16]. The sensitivities in these ranges will not be studied in the framework of this work.

1.2 Evaluation of Nuclear Data

Resonance parameters and cross sections are part of so-called nuclear data that are generally compiled into libraries. The cross sections used in stochastic and deterministic codes are retrieved from these libraries: For example, the XDR files used by TRIPOLI-4® _{are generated with these} libraries. Depending on the user purpose, these cross sections could be defined by groups of energy, or by pointwise range. Regularly, these libraries are updated to improve their accuracy and consistency.

1.2.1 Evaluated Nuclear Data File

Different international libraries are used across the world as nuclear data evaluation references, like: the Joint Evaluated Fission and Fusion File (JEFF, the European library), Evaluated Nuclear Data File (ENDF, the American library) and Japanese Evaluated Nuclear Data Library (JENDL). The different actors behind these libraries have decided to edit then in an international format, ENDF-6 [14], based on the ENDF library format, and currently in its 6th version.

These libraries are regularly upgraded and improved, by direct experiments and integral experiments assimilations. P. Tamagno [15] have illustrated in Figure 1.3 the complexity of upgrading nuclear libraries:

It should be noticed that evaluated data are not standards, which are constructed from specific microscopic experiments, designed to study a particular reaction at an exclusive energy range. Such libraries are based on theoretical model fitted on experimental information, and the library consistency is improved using the feedback of integral experiment. Indeed, microscopic experimental information have a restricted domain of energy, and sometimes are scarce [17].

24 | Evaluation of cross sections

A library contains a file for each used material at a reference temperature of 0 K (the Doppler broadening is generally done during the processing). Each file is composed by subdivision called MF, and each MF is split by reaction, one reaction being associated with one subdivision called MT.

The following tables summarize some of the main MF and MT:

Table 1.1 - Extract of various MF [10]

MF Description

1 General information 2 Resonance parameter data 3 Reaction cross sections

4 Angular distributions for emitted particles 5 Energy distributions for emitted particles … …

32 Data covariance for resonance parameters

Table 1.2- Extract of various MT [10]

MT Description

1 Total cross sections (incident neutrons only) 2 Elastic Scattering Cross Section

16 MT for emission of 2 neutrons (in addition to the residual nucleus) 18 Fission

50-90 Neutron emission for each nucleus states (50=ground state) 102 Radiative Capture

Data stored in each MT can be pointwise parameters (cross sections in MF3), but also characteristics of distribution law (Legendre polynomial coefficients in MF4). Some possess flags to indicate the law/representation used to evaluate the parameters (LRF flag in MF2, LF flag in MF5). As pointed in the Table 1.1, data about uncertainties and correlations could also be found for a MF/MT couple in an ENDF file. Even if these libraries are nuclear data references, they still need improvements to counter 5 major aspects, each of them being pointed and exemplified in [15]:

• Consistency of nuclear files (discontinuity between URR and Continuum model)

• Consistency of nuclear data modelling (the same isotope could have different parameters in different reaction)

• Lack of variance-covariance data

• Predictability of not-measured nuclear data (experiments are hardly practical, expensive and need alternatives)

• Capability to reproduce experimental data

Evaluation of cross sections | 25

1.2.2 CONRAD

In the framework of the JEFF project [18], the LEPh develops the code CONRAD since 2005 [19]. This is an object-oriented C++ code, which can work with various formalisms in Resolved Resonances Range (RRR) (R-Matrix approximation, Reich-Moore, Multi-Level Breit-Wigner), Unresolved Resonances Range (URR) (Average R-Matrix) and continuum (Optical model). It is particularly used during the integral experiments assimilation process, as its functionalities cover theoretical calculation of cross section, adjustment of models on experimental data, Doppler broadening ... The structure of CONRAD consists in a dynamic library linked to an executable. The library is split into different classes, which cover analysis tools, theoretical models, experimental aspects, interfaces...

This project will firstly use the available ENDF parser as a basis to reconstruct ENDF files. Subsequently, CONRAD abilities to reconstruct cross sections and calculate resonance parameters derivatives will be needed for the development of the IFP method. Indeed, CONRAD possesses various classes to:

• Reconstruct cross sections and calculates their resonance parameters derivatives, from given resonance parameters and an energetic grid. These classes will construct the collision matrix for the formalism required and calculate the cross section at each point of the given energetic grid.

• Use the free gas model to broaden the cross section at the desired temperature and calculate the new derivatives.

• Determine the optimal energetic grid to represent the cross sections.

1.3 Summary

Exact perturbations in Monte Carlo simulation | 27

2 Exact perturbations in Monte Carlo simulation

The previous chapter have highlighted two major facts: cross sections can be reconstructed from resonance parameters, and, moreover, cross sections sensitivities to resonance parameters can be calculated with the help of CONRAD. This chapter details how reactivity sensitivities to cross sections are computed in this work.

2.1 Monte Carlo simulations

Monte Carlo simulations are often used to calculate reference integral parameters, during a reactor design study for example. This part will briefly present the principle of the Monte Carlo method, and the reference code used at the CEA: TRIPOLI-4®_.

2.1.1 Principle

Monte Carlo simulations are based on a stochastic approach. In the framework of nuclear engineering, they are used to obtain a stochastic solution to a deterministic problem: solving of the neutron transport equation [20]. This method gives a mean value that should be equal to the expected value of the determinist approach.

More precisely, the code follows a particle life: for example, in a criticality calculation, the code will generate neutrons with pre-set initial direction, energy and location. It will then use probability distributions to determine the next position and interaction of the neutron, and so identify whether the neutron is terminated (absorption, leakage, out of considered energy range), or will continue (scattering).

In criticality calculation, such process will be iterated on a large number of neutrons to estimate the reactivity of the studied system, using the effective neutron multiplication factor keff, usually defined as the ratio of a neutron population at generation n by the population at generation n-1. This definition could be slightly different in TRIPOLI-4®_{, depending on the estimator used (see} Section 2.1.2).

The current benefit of stochastic approach is that approximations used in deterministic solving could be overpassed. On the other hand, the computing effort needed is a major drawback.

2.1.2 TRIPOLI-4®

TRIPOLI-4® _{is the reference Monte Carlo code for the three major French industrial stakeholders} (CEA, AREVA and Electricité de France (EDF)). Developed by the Service for Reactor Studies and Applied Mathematics (SERMA), it is the fourth generation of the CEA Monte Carlo code (since the mid-60) [21]. It manages large parallelization to counteract the disadvantage of Monte Carlo simulations. Diverse kinds of particles could be tracked: neutrons of course (20 MeV till 10 −5 eV), but also photons, electrons, positrons (at different energy range). It offers source-fixed as well as criticality simulation modes, and most of the common variance reduction technics, as the Russian roulette.

28 | Exact perturbations in Monte Carlo simulation

decreased. For example, for a capture, instead of killing the neutrons, its weight will be reduced. When the weight is too low, the neutron history is terminated. With an appropriated correction, the effective neutron multiplication factor can be found. As the number of interactions sampled is larger than for a standard Monte Carlo simulation, the variance of the result is improved.

Regarding the estimators of the effective neutron multiplication factor, TRIPOLI-4® _{uses three} different estimators [22]:

• kSTEP, which considers the neutrons the code produces at each generation

• kCOLL, which considers the average neutrons the code should produce, based on the collisions sampled.

• kTRACK, which considers the average neutrons the code should produce, based on the neutrons tracks.

Sensitivities studies are also available in TRIPOLI-4®_{, using correlated sampling method or} Taylor series development (both giving approximated sensitivities). Recent version possesses sensitivities calculation implemented by G.Truchet [3], both based on the Iterated Fission Probability (IFP): the Collisions-based Exact Perturbations (CEP) and the Black Body Exact Perturbations (BBEP) methods.

2.2 Adjoint Flux

The next part, dedicated to the calculus of the reactivity sensitives, requires the introduction of a specific tool frequently used in deterministic codes, and recently implemented in several Monte Carlo code [2] [23]: the adjoint flux.

2.2.1 Definition

Considering two real functions 𝜙 and 𝜓:

{𝜙 ∶ (𝑥⃗, Ω, 𝐸, 𝑡) → ℝ

𝜓 ∶ (𝑥⃗, Ω, 𝐸, 𝑡) → ℝ (2.1)

The scalar product of these two functions is:

〈𝜓, 𝜙〉 = ∫ ∫ ∫ 𝜓(𝑥⃗, Ω, 𝐸)𝜙(𝑥⃗, Ω, 𝐸)𝑑Ω𝑑𝑉𝑑𝐸 𝜋 𝑅 ∞ 0 (2.2)

Adjoint operator ℒ† _{from a given linear operator ℒ is defined by the following property, for all 𝜙} and 𝜓:

〈𝜓, 𝓛𝜙〉 = 〈𝜙, 𝓛†_{𝜓〉 (2.3)}

Exact perturbations in Monte Carlo simulation | 29

𝓗𝜙 = 0 (2.4)

The 𝓗 operator could be split into two: one operator referring to the neutrons production 𝓟, and a second to the neutrons disappearance 𝓚, so that for a critical and stationary reactor:

𝓗𝜙 = 𝓟𝜙 − 𝓚𝜙 = 0 (2.5)

If a reactor is not critical, it could still be described by the previous equation, by balancing 𝓟 with a factor 𝜆: For an overcritical reactor, 𝜆 will be used to mathematically reduce the term of neutrons production, and conversely. 𝜆 is not more than the opposite of the reactor multiplication factor.

Using Equation (2.4):

𝓚𝜙 = 𝜆𝓟𝜙 (2.6)

The adjoint operators of 𝓚 and 𝓟 exist [3], and the function 𝜙†_{, called adjoint flux, verify} with 𝜙:

{〈𝜙

†_{, 𝓟𝜙〉 = 〈𝜙, 𝓟}†_𝜙†_〉

〈𝜙†_{, 𝓚𝜙〉 = 〈𝜙, 𝓚}†_𝜙†_〉 (2.7)

It could easily be proved that 𝜆† _{= 𝜆 as the flux and its adjoint could not physically be null.} 2.2.2 Iterated Fission Probability

Physically, the adjoint flux could be seen as a characterization of a neutron importance: Depending on its initial position and its energy (thermal or fast), a neutron will have a higher probability to father a fission. The importance could then be quantified by evaluating the asymptotic average descendants of a neutron: The higher the importance of a neutron will be, the more neutrons it will generate in the long term. This importance is called Iterated Fission Probability (IFP). Mathematically, a proportional relationship could be established between the IFP and the adjoint flux [1].

For a neutron, the number of descendants N it will father after L generations (so-called IFP cycle) could be estimated by [3]:

𝑁𝐿_{= ∏ 𝑘𝑔} 𝐿

𝑔=1

(2.8)

𝑘𝑔 being the keff of the generation g. Of course, for an over-critical reactor, such a number will tend to infinity. The importance I after L generations is then defined as [3]:

𝜙†_{∝ 𝐼}𝐿_{= ∏} 𝑘𝑔 𝑘𝑔→∞ 𝐿

𝑔=1

30 | Exact perturbations in Monte Carlo simulation

This definition shows the need to previously carry out an accurate and well converged Monte-Carlo simulation before computing the adjoint flux, to obtain an estimation of 𝑘𝑔→∞ . Indeed, the IFP method implemented by G. Truchet distorts the calculation of the effective neutron multiplication factor in TRIPOLI-4®_{. Other methods exist to compute the adjoint flux [3], but none} of them are implemented in TRIPOLI-4®_.

2.3 Reactivity sensitivities

This part will focus on the two methods available to calculate reactivity sensitivities, how G. Truchet implemented the Collision based Exact Perturbation (CEP) method in TRIPOLI-4®_{, and how} sensitivities could be computed from this method.

2.3.1 Eigenvalue difference method

Considering a given parameter Χ whose influence on a given system’s reactivity should be studied. Here, this parameter could be, for instance, a cross section, a nuclide density or a resonance parameter. The most common way to calculate reactivity sensitivities to this parameter is to use “direct” calculation. Such a method consists in doing two Monte Carlo simulations: One in a reference state, a second were the parameter is slightly modified. The sensitivity (expressed in pcm/%) can then be easily calculated with a finite difference, using a first order Taylor series:

𝜕𝜌 𝜕Χi≈ 𝛿𝜌 𝛿Χ𝑖= 𝜌𝑝𝑒𝑟𝑡− 𝜌𝑟𝑒𝑓 𝜒_{𝑖 𝑝𝑒𝑟𝑡}− Χ_{𝑖 𝑟𝑒𝑓}= 1 Χ_{𝑖 𝑝𝑒𝑟𝑡}− Χ_{𝑖 𝑟𝑒𝑓}⋅ ( 1 𝑘𝑟𝑒𝑓− 1 𝑘𝑝𝑒𝑟𝑡) (2.10)

Neglecting correlation between simulations, the sensitivity variance is given by:

𝜎² (𝜕𝜌 𝜕Χi ) = ( 1 Χ_{𝑖 𝑝𝑒𝑟𝑡}− Χ_{𝑖𝑟𝑒𝑓}) 2 ⋅ [𝜎² ( 1 𝑘𝑟𝑒𝑓 ) + 𝜎² ( 1 𝑘𝑝𝑒𝑟𝑡 )] (2.11)

This expression could be developed, for 𝜎[𝑘𝑖] ≪ 𝑘𝑖, with a first order Taylor series, as: 𝜎² (𝜕𝜌 𝜕Χi ) ≈ ( 1 Χ_{𝑖𝑝𝑒𝑟𝑡}− Χ_{𝑖 𝑟𝑒𝑓}) 2 ⋅ [𝜎²(𝑘𝑟𝑒𝑓) 𝑘𝑟𝑒𝑓4 +𝜎²(𝑘𝑝𝑒𝑟𝑡) 𝑘𝑝𝑒𝑟𝑡4 ] (2.12)

This approximation of the sensitivity variance is largely sufficient for the analysis conducted in this thesis: Indeed, tests conducted by G. Truchet [3] have shown that for 𝑘𝑟𝑒𝑓 one, a 5000 pcm reactivity difference would only induce a bias of 1% in the variance estimated. These conditions are always verified for reactivity difference encountered chapter five.

Exact perturbations in Monte Carlo simulation | 31

2.3.2 Exact and first order Perturbation Theory

Considering two different reactors states, named ref and pert, it should be kept in mind that: Δ𝜌 = 𝜌𝑝𝑒𝑟𝑡− 𝜌𝑟𝑒𝑓=

1 𝑘𝑝𝑒𝑟𝑡−

1

𝑘𝑟𝑒𝑓 = 𝜆𝑝𝑒𝑟𝑡− 𝜆𝑟𝑒𝑓 (2.13)

With 𝜆𝑖 solution of the Equation (2.6) for a state 𝑖. Using the Equations (2.6), for a reference (1) and a perturbed state (2):

{(𝜆1𝓟1− 𝓚1)Φ1= 0 (𝜆2𝓟2− 𝓚2)Φ2= 0

(2.14)

By substituting the first with the second, and extracting a factor (𝜆1− 𝜆2) [20]:

(𝜆1− 𝜆2)𝓟1𝜙1 = (𝜆2𝓟2− 𝓚2)(𝜙2− 𝜙1) + [(𝜆2( 𝓟2− 𝓟1) − (𝓚2− 𝓚1)]𝜙1 (2.15) Now, the scalar product with 𝜙₂† _{give us:}

(𝜆1− 𝜆2)〈 𝜙2 †

, 𝓟1𝜙1〉 = 〈𝜙₂†, (𝜆2𝓟2− 𝓚2)(𝜙2− 𝜙1)〉 + 〈𝜙2 †

, [(𝜆2( 𝓟2− 𝓟1) − (𝓚2− 𝓚1)]𝜙1〉 (2.16) As seen in (2.3), and applying (2.6)

〈𝜙₂†, (𝜆2𝓟2− 𝓚2)(𝜙2− 𝜙1)〉 = 〈(𝜆2𝓟2 †_{− 𝓚} 2 †_)𝜙 2 †_{, (𝜙} 2− 𝜙1)〉 = 0 (2.17)

Equation (2.16) could then be condensed in:

Δ𝜌 = (𝜆1− 𝜆2) =〈𝜙2 †

, (𝜆2Δ𝓟 − Δ𝓚)𝜙1〉 〈 𝜙₂†, 𝓟1𝜙1〉

(2.18) Considering small perturbations [26], the perturbed adjoint flux and production operator could be defined as: {𝜙2 † = 𝜙₁†+ 𝛿𝜙₁† 𝓟2= 𝓟1+ 𝛿𝓟1 (2.19)

By only keeping the first order terms, 𝛿𝜌 is immediately given by:

𝛿𝜌 = 〈𝜙1 †

, (𝜆1𝛿𝓟 − 𝛿𝓚)𝜙1〉 〈 𝜙₁†, 𝓟1𝜙1〉

(2.20) Many advantages come from (2.20): Flux in the perturbed state are not required to obtain the reactivity variation from any kind of perturbation (concentration, cross section, geometry …). Flux should also be computed only once for diverse types of perturbation.

2.3.3 Collision based Exact Perturbation (CEP)

32 | Exact perturbations in Monte Carlo simulation 𝓚𝜙1(𝑟,⃗⃗⃗ 𝐸, Ω) = 𝑑𝑖𝑣(Ω ⋅ 𝜙1(𝑟,⃗⃗⃗ 𝐸, Ω)) + Σ𝑡(𝑟,⃗⃗⃗ 𝐸) ⋅ 𝜙1(𝑟,⃗⃗⃗ 𝐸, Ω) − ∫ Σ𝑠(𝑟⃗ 𝐸′_,Ω′ , 𝐸′_{→ 𝐸, Ω}′_{→ Ω) ⋅ 𝜙} 1(𝑟,⃗⃗⃗ 𝐸′, Ω′)𝑑𝐸′𝑑²Ω′ _(2.21)

And the production operator 𝓟, with 𝜈 the number of neutrons: 𝓟𝜙1(𝑟,⃗⃗⃗ 𝐸, Ω) =

1

4𝜋∫_𝐸′_,Ω′𝜈(𝑟,⃗⃗⃗ 𝐸′) ⋅ Σ𝑓(𝑟,⃗⃗⃗ 𝐸′) ⋅ 𝜒𝑓(𝑟⃗

, 𝐸′_{→ 𝐸) ⋅ 𝜙}

1(𝑟,⃗⃗⃗ 𝐸′, Ω′)𝑑𝐸′𝑑²Ω′ (2.22) The variation of these operators could then be expressed as:

Δ𝓚𝜙1(𝑟,⃗⃗⃗ 𝐸, Ω) = ΔΣ𝑡(𝑟,⃗⃗⃗ 𝐸) ⋅ 𝜙1(𝑟,⃗⃗⃗ 𝐸, Ω) − ∫ ΔΣ𝑠(𝑟⃗ 𝐸′_,Ω′ , 𝐸′_{→ 𝐸, Ω}′_{→ Ω) ⋅ 𝜙} 1(𝑟,⃗⃗⃗ 𝐸′, Ω′)𝑑𝐸′𝑑²Ω′ (2.23) Δ𝓟𝜙1(𝑟,⃗⃗⃗ 𝐸, Ω) = 1 4𝜋∫𝐸′_,Ω′Δ[𝜈(𝑟,⃗⃗⃗ 𝐸′) ⋅ Σ𝑓(𝑟,⃗⃗⃗ 𝐸′) ⋅ 𝜒𝑓(𝑟⃗, 𝐸 ′_{→ 𝐸)] ⋅ 𝜙1(𝑟,⃗⃗⃗ 𝐸}′_{, Ω}′_{)𝑑𝐸′𝑑²Ω′ (2.24)} After multiplying by 𝜙2 †_(𝑟,

⃗⃗⃗ 𝐸, Ω) and integrating on each side over all 𝑟⃗, 𝐸 and Ω Equations (2.23) and (2.24): 〈𝜙₂†, Δ𝓚𝜙1〉 = 〈𝜙₂†, ΔΣ𝑡𝜙1〉 − 1 4𝜋∫ 𝜙2 †_{(𝑟,⃗⃗⃗ 𝐸}′_{, Ω}′₎ 𝑟⃗,𝐸′_,Ω′ ∫ ΔΣ𝑠(𝑟⃗ 𝐸,Ω , 𝐸 → 𝐸′_{) ⋅ 𝜙} 1(𝑟,⃗⃗⃗ 𝐸, Ω)𝑑𝐸𝑑²Ω𝑑𝐸′𝑑²Ω′𝑑3𝑟 ⃗⃗⃗ (2.25) 〈𝜙₂†, Δ𝓟𝜙1〉 = 1 4𝜋∫ 𝜙2 † (𝑟,⃗⃗⃗ 𝐸′, Ω′) 𝑟⃗,𝐸′_,Ω′ ∫ Δ[𝜈(𝑟,⃗⃗⃗ 𝐸) ⋅ Σ𝑓(𝑟,⃗⃗⃗ 𝐸) ⋅ 𝜒𝑓(𝑟⃗ 𝐸,Ω , 𝐸 → 𝐸′_{)] ⋅ 𝜙} 1(𝑟,⃗⃗⃗ 𝐸, Ω)𝑑𝐸𝑑²Ω𝑑𝐸′𝑑²Ω′𝑑3𝑟 ⃗⃗⃗ (2.26)

Because of the scattering and the fission, it appears immediately that 〈𝜙2 †

, (𝜆2Δ𝓟 − Δ𝓚)𝜙1〉 could not be decomposed as a sum of the scalar products 〈𝜙2

† , ΔΣ𝑐𝜙1〉, 〈𝜙2 † , ΔΣ𝑓𝜙1〉, 〈𝜙2 † , ΔΣ𝑠𝜙1〉. Such an expression will have made the perturbation computation simpler. Nevertheless, the Collision based Exact Perturbation [3] method offers a way to overpass this difficulty, and to compute 𝓚 and 𝓟 using neutron collisions.

This method consists in storing various parameters for each neutron collision 𝑖, as: the weight 𝑤𝑖 of the incident neutron, its energy, position, direction. The cross sections concerned by the collision are also stored. These parameters are stored in so-called “Stock files”. In a second run, flux is calculated with IFP for each collision, and the operators’ variations Δ𝓟 and Δ𝓚 are calculated. This way, the total reactivity perturbation could be found by summing the weighted contribution of all collisions. The flux estimator at a collision 𝑖 is expressed with 𝑤𝑖 the weight of the incident neutron and Σ𝑡,𝑖 the total cross section at the collision I by [3]:

𝜙𝑖= 𝑤𝑖

Σ𝑡,𝑖 (2.27)

From now on, the adjoin flux estimator in the perturbed state is the importance, named IL_. Equations (2.25) and (2.26) show that the product 〈𝐼𝐿,(𝜆2Δ𝓟−Δ𝓚)𝜙〉

〈𝐼𝐿_,𝓟𝜙〉 can be decomposed into a sum of each reaction type contribution:

〈𝐼𝐿_{, (𝜆2Δ𝓟 − Δ𝓚)𝜙〉}

Exact perturbations in Monte Carlo simulation | 33

With the CEP method, each of these Δ𝜌 could be computed by summing the contribution of each collision:

Δ𝜌𝑐 = ∑ δ𝜌𝑐,𝑖 𝑖

(2.29) Taking, for instance, a nuclide n, one collision 𝑖, at the position 𝑟 ⃗⃗⃗, with an incident neutron at the energy 𝐸 and a direction Ω. The contribution 𝛿𝜌𝑐,𝑖 could be expressed as:

𝛿𝜌𝑛,𝑐,𝑖 = −I𝑖𝐿 (𝑟,⃗⃗⃗ 𝐸, Ω) ⋅ ΔΣ𝑛,𝑐(𝑟,⃗⃗⃗ 𝐸, Ω) ⋅ 𝜙𝑖(𝑟,⃗⃗⃗ 𝐸, Ω) (2.30) For a simple case, where the incident neutron with an energy and direction (𝐸, Ω) can only scatter to (𝐸′, Ω′), the contribution 𝛿𝜌𝑛,𝑠,𝑖 could be expressed as

𝛿𝜌𝑛,𝑠,𝑖 = (I𝑖𝐿(𝑟,⃗⃗⃗ 𝐸′, Ω′) − I𝑖𝐿(𝑟,⃗⃗⃗ 𝐸, Ω) ) ⋅ 𝜙𝑖(𝑟,⃗⃗⃗ 𝐸, Ω) ⋅ ΔΣ𝑛,𝑠(𝑟⃗, 𝐸 → 𝐸′, Ω → Ω′) (2.31)

With the example depicted in Figure 2.1, it could immediately be seen that for the collision 𝑖, the importance of the neutron in the state (𝑟,⃗⃗⃗ 𝐸, Ω) is needed, as well as the importance of the neutron in all the states (𝑟,⃗⃗⃗ 𝐸′_{, Ω}′_{) possible. The same reasoning could be applied to fission reactions. In} practice, collisions are stored after a direct calculation in the reference state. Thereafter, an adjoint flux calculation is done in the perturbed state for each possible subsequent reaction, with an new Monte Carlo simulation. Summing the contribution of all collisions, the reactivity variation can be computed using the perturbation theory, as all the parameters needed are available in the Stock Files.

Collision files are firstly generated in TRIPOLI-4® _{and all the perturbation calculations are} made in the post-processing code JIMMY2. To reduce computing time and the resources used, P. Tamagno have recently implemented an improved version of JIMMY2, that converts Stock Files into lighter Files, as it does not contain the collisions cross sections and concentrations. These are retrieved from libraries and inputs files during the post-processing.

34 | Exact perturbations in Monte Carlo simulation

2.3.4 Sensitivities in JIMMY2

Reactivity sensitivities are, in JIMMY2, expressed in pcm per percent. For a reaction r, a nuclide n:

𝑆𝜌,𝜎𝑛,𝑟= ∫ 𝜎𝑛,𝑟(𝑟,⃗⃗⃗ 𝐸) ⋅ 𝜕𝜌(𝑟,⃗⃗⃗ 𝐸, Ω) 𝜕𝜎𝑛,𝑟(𝑟,⃗⃗⃗ 𝐸)𝑑𝐸𝑑²𝛺𝑑 3_{𝑟 ⃗⃗⃗} 𝑟, ⃗⃗⃗⃗𝐸,Ω ≈ ∫ 𝜎𝑛,𝑟(𝑟,⃗⃗⃗ 𝐸) ⋅ 𝛿𝜌(𝑟,⃗⃗⃗ 𝐸, Ω) 𝛿𝜎𝑛,𝑟(𝑟,⃗⃗⃗ 𝐸) 𝑟, ⃗⃗⃗⃗𝐸,Ω 𝑑𝐸𝑑²𝛺𝑑3_{𝑟 ⃗⃗⃗ (2.32)}

Applying (2.32) on the capture section of a nuclide n: 𝑆𝜌,𝜎_𝑛,_𝑐≈ ∫ 𝜎𝑛,𝑐(𝑟,⃗⃗⃗ 𝐸) ⋅ 𝛿𝜌(𝑟,⃗⃗⃗ 𝐸, Ω) 𝛿𝜎𝑛,𝑟(𝑟,⃗⃗⃗ 𝐸) 𝑟, ⃗⃗⃗⃗𝐸,Ω 𝑑𝐸𝑑²𝛺𝑑3𝑟 ⃗⃗⃗ = ∫ 𝜎𝑛,𝑐(𝑟,⃗⃗⃗ 𝐸) ⋅ ∑ 𝛿𝜌𝑖 𝑛,𝑐,𝑖(𝑟,⃗⃗⃗ 𝐸, Ω) 𝛿𝜎𝑛,𝑐(𝑟,⃗⃗⃗ 𝐸) 𝑟, ⃗⃗⃗⃗𝐸,Ω 𝑑𝐸𝑑²𝛺𝑑3𝑟 ⃗⃗⃗ _(2.33) Indeed, only the capture contribution to reactivity variation will influence sensitivity of reactivity to capture cross section. Thereby, decompositions seen in Equations (2.28) and (2.29) are applied. Using Equation (2.30) and the approximation for small perturbations from Equation (2.20): ∫ ∑ 𝜎𝑛,𝑐(𝑟 ⃗⃗⃗, 𝐸) 𝛿𝜎𝑛,𝑐(𝑟⃗, 𝐸)⋅ I𝑖𝐿(𝑟 ⃗⃗⃗, 𝐸, Ω) ⋅ 𝛿Σ𝑛,c(𝑟,⃗⃗⃗ 𝐸) ⋅ 𝛿𝜙𝑖(𝑟 ⃗⃗⃗, 𝐸, Ω) 〈IL_{, 𝓟𝜙〉} 𝑖 𝑟 ⃗⃗⃗⃗,𝐸,Ω 𝑑𝐸𝑑²𝛺𝑑3𝑟 ⃗⃗⃗ = ∫ ∑ 𝜎𝑛,𝑐(𝑟 ⃗⃗⃗, 𝐸) 𝛿𝜎𝑛,𝑐(𝑟 ⃗⃗⃗, 𝐸) ⋅ I𝑖 𝐿_{(𝑟 ⃗⃗⃗, 𝐸, Ω) ⋅ 𝑁𝑛⋅ 𝛿𝜎𝑛,𝑐(𝑟 ⃗⃗⃗, 𝐸) ⋅ 𝛿𝜙(𝑟 ⃗⃗⃗, 𝐸, Ω)} 〈IL_{, 𝓟𝜙〉} 𝑖 𝑟 ⃗⃗⃗⃗,𝐸,Ω 𝑑𝐸𝑑²𝛺𝑑3_{𝑟 ⃗⃗⃗} = 〈I L_{, 𝑁} 𝑛⋅ 𝜎𝑛,𝑐⋅ 𝜙〉 〈IL_{, 𝓟𝜙〉} (2.34)

The same process can be achieved on the scattering and fission reaction. As the approximation of the sensitivity is done using Equation (2.20), the adjoint calculation can be done on the reference state only, as depicted in Figure 2.2. It should be recalled that the adjoint flux value at each collision rely on the reaction studied (see Section 2.3.3).

Exact perturbations in Monte Carlo simulation | 35

2.4 Summary

Methodology | 37

3 Methodology

The points discussed previously have shown that cross sections are dependent on resonance parameters, and that sensitivity of reactivity to cross sections can be computed with at least two methods. This chapter will discuss how to extend these methods to compute sensitivities of reactivity to parameters resonances.

3.1 Sensitivities to resonance parameters

From the Section 2.3, three methods have been defined to achieve sensitivities to resonance parameters computation:

1. As introduced in 2.1.1, the easiest way to obtain these sensitivities is to proceed by finite difference. This is achieved by using two evaluation files: one unchanged, to simulate the reference state with a “direct” Monte Carlo computation; and one with a modified resonance parameter, to simulate the perturbated state. Sensitivities to resonance parameters are then calculated by:

𝜕𝜌 𝜕Γi ≈𝛿𝜌 𝛿Γ𝑖 = 𝜌𝑝𝑒𝑟𝑡− 𝜌𝑟𝑒𝑓 Γ_{𝑖𝑝𝑒𝑟𝑡}− Γ_{𝑖𝑟𝑒𝑓} = 1 Γ_{𝑖𝑝𝑒𝑟𝑡}− Γ_{𝑖𝑟𝑒𝑓}⋅ ( 1 𝑘𝑟𝑒𝑓 − 1 𝑘𝑝𝑒𝑟𝑡 ) (3.1)

2. The second way is to calculate exact perturbation using the original IFP method described in Section 2.3.2. The direct computation is done using the unchanged evaluation file, and the adjoint running with the modified one. The major difference with the “direct” computation is that convergence is now made on the reactivity variation, and not the reactivity itself: 𝛿𝜌 is directly obtained after the simulation.

This method can improve the computing time, but the process is laborious: completely new computation and TRIPOLI-4® _{libraries XDR files are needed for each parameter studied. It} will be used only to test the original IFP method running.

3. The last possibility is the one developed in this project: Calculate sensitivity to resonance parameters directly inside the IFP method. For this purpose, CONRAD abilities to calculate parameters resonances derivatives are implemented in the post-processing JIMMY2. This way, there are no more needs to modify library files, and the collision files should not be recomputed for different nuclei or parameters.

3.2 Approach

To ensure the reliability of the sensitivities calculated with the new IFP method, reference sensitivities must be found using direct calculation. For this purpose, TRIPOLI-4® _{libraries must be} produced with possibly modified parameters. These files will also be used to test the original IFP method: If equivalent sensitivities are found, the second method can be safely used with CONRAD sensitivities to develop the third method. Following steps will be executed:

• Improvement and test of the current ENDF parser in CONRAD, to rebuild complete ENDF files with modified MF2.