Go To: "SMURE Class list"

Report: LPSC 17002

MURE v.3 : SMURE,

Serpent-MCNP Utility for Reactor Evolution

User Guide - 2nd Edition

pict

January 2024

Main Contributors :

See also the FAQ

Contents

1 Introduction

2 What’s new in MURE

3 From MURE to SMURE : Choosing the Monte-Carlo Transport Code

4 Geometry Definition

5 Materials, Sources, Tallies, ...

6 Nuclei Tree

7 Evolution in MURE

8 Looking at results from the evolution

9 Thermal-hydraulics/neutronics coupling

A Basic of C++ to understand MURE

B Node tree simplification

C Stochastic volume calculation

D Back-Stage processes: everything you always wanted to know about MURE* (*But were afraid to ask)

Bibliography

Chapter 1
Introduction

In the following, we refer either to MURE (MCNP Utility for Reactor Evolution), MURE ≥ 2 or SMURE (Serpent/MCNP Utility for Reactor Evolution)1 ; they represent the same package (MURE is the version 1 of SMURE-MURE ≥2). The main aim of the MURE[1, 2]/SMURE package is to perform nuclear reactor time-evolution using the widely-used particle transport code MCNP[3] (a Monte Carlo code which is mostly written in FORTRAN) or Serpent2[4, 5]. Many depletion codes exist for determining time-dependent fuel composition and reaction rates. These codes are either based on solving Boltzman equation using deterministic methods or based on Monte-Carlo method for neutron transport. Among them, one has to cite MCNPX/CINDER 90[6], MONTEBURN[7], KENO/ORIGEN[8], MOCUP[9], MCB[10], VESTA/MORET[11, 12], TRIPOLI-4D[13], Serpent[4], ... which provide neutron transport and depletion capabilities. However, the way to control (or interact with) the evolution are either limited to specific procedure and/or difficult to implement.

In (S)MURE, due to the Object-oriented programming, any user can define his own way to interact with evolution. From an academic point of view, it is also good to have lots of M-C evolution codes to compare and benchmark them to understand physics approximations of each one. Moreover, SMURE provides a simple graphical interface to visualize the results. It also provides a way to couple the neutronics (with or without fuel burn-up) and thermohydraulics using either an open source simple code developed in SMURE (BATH, Basic Approach of Thermal Hydraulics) or a sub-channel 3D code, COBRA-EN[14, 15]. But SMURE can also be used just as an interface to MCNP or Serpent to build geometries (e.g. for neutronics experiments simulation).

SMURE is based on C++ objects allowing a great flexibility in the use2. There are 4 main parts in this library:

  1. Definition of the geometry, materials, neutron source, tallies, ...

  2. Construction of the nuclear tree, the network of links between neighbouring nuclei via radioactive decays and nuclear reactions.

  3. Evolution of some materials, by solving the corresponding Bateman’s equations.

  4. Thermal-hydraulics: it couples neutronics, thermal-hydraulics and, if needed, fuel evolution.

The use of this package requires the following installation:

  1. a C++ compiler (mandatory). SMURE is developed using gcc (for the more recent implementation the use of c++14 standard is required).

  2. Serpent2, MCNP or MCNPX (mandatory, available at the NEA DataBank & RSICC). These codes will be referred as “MC” in the following.

    1. Nuclear data files process for MCNP (ACE format) and/or NJOY[16] to process ENDF files into the ACE format (Nuclear data in ACE format are mandatory).

    2. Serpent 1 can not be used with MURE because it does not support unions.

  3. COBRA-EN for the thermal-hydraulic part (available the atNEA DataBank & RSICC)

  4. The ROOT graphical tools developed at CERN (http://root.cern.ch) is necessary if user wants to use the post treatment GUI tools. In the GUI, radiotoxicity of fuel/waste can be calculated and plotted if the LAPACK library is installed (Linear Algebra PACKage, see http://www.netlib.org/lapack/index.html or standard LINUX repositories (Redhat, Fedora, Ubuntu, ...) for pre-installed versions). If you are using ROOT-5 or before, you should modify the beginning of the Makefile of MureGui.

At present, SMURE has only been compiled and tested on LINUX/UNIX platforms.

1.1 Installation

TO BE NOTICED: for some evident reasons5, upgrading for MURE version 1 to MURE version 2 breaks some backward compatibility. Nevertheless, for “old” MURE user changes will be minors ; we have try to keep the essential of MURE, but given more flexibility. A chapter is devoted to the migration of old MURE users.

1.1.1 Compilation

Remak:

1.1.2 Remarks on MureGui compilation

The compilation is enable if you invoke “build.sh” script with “--enable-gui” or if you are using directly the “cmake” command, you should use the -DUSE-GUI=ON” flag.. Using MureGui without argument gives a short description of the code used (see also § 8.3). You need a modern version of ROOT (ROOT 5 or later). Optional libraries for radiotoxicity (LAPACK & BLAS) are also needed (in general, you need to install the “lapack-devel” on most of Linux distributions) ; installation of ROOT and these libraries are needed BEFORE running the “build.sh” script or cmake (or you have to re-run it).

1.1.3 Remarks on MCNP/makxsf compilation

In MURE, the access to ACE MCNP cross-sections is done for different purposes (find if cross-sections exist, find the total cross-sections in order to write in MC input file only “significant nucleus” (see 5.1.6), use of “multi-group flux” evolution method (see item 5 of section 7.3), ...). The use of binary ACE files improves the reading time (and also the disk space necessary to store the cross-sections). In C/C++, reading/writing in unformatted files (binary) is done by reading records of 1 byte long. But some FORTRAN compilers such as IFORT of Intel read/write in unformatted file with 4 bytes records. Thus if you want to use such a compiler for MCNP/makxsf you have to force the record length in unformatted files to be 1 byte. In the Intel IFORT compiler this is done via the –assume byterecl flag of ifort command line (see the System.gcf file of the config dir of MCNP distribution where System is either Linux, AIX, …).

1.1.4 Remarks on MURE with Serpent2

At the present time, Serpent2 seems only to handle ASCII file for ACE cross-section. Thus one has to use this ASCII data if one use Serpent. Conversion from the MCNP xsdir to Serpent xsdata is done either with the code provided with Serpent, either directly in MURE if one use the MURE:SetAutoXSDIR() method.

1.1.5 Building files for evolution

In this section, it is explained how to build the 2 necessary files BaseSummary.dat and AvailableReactionChart.dat (these files are not provided in the distribution). One supposes that a user has a standard MCNP6 ENDF B6 base provided with MCNP6. The provided xsdir file and the nuclear data are located in /path_2_ace_files directory (Verify that in the xsdir the absolute path to data files is present).

1    1    0 ENDFB  6.1  293.62  .24c  STD    1001.24c 0.999170 /path_2_ace_files/la150n2 0 2 1 10106 4096 512 2.5301E-08 
...
1    1    0 ENDFB  6.1  293.62  .66c  STD    1001.66c 0.999170 /path_2_ace_files/endf66a2 0 2 1 10128 4096 512 2.5301E-08 
...
92   235  0 ENDFB  6.1  293.62  .66c  STD  92235.66c 233.025000 /path_2_ace_files/endf66c2 0 2 6899 722105 4096 512 2.5301E-08 ptable 
...

1.1.6 Running some examples

The directory MURE/examples contains commented examples on how to use MURE. When this is possible, the same example exist for both MCNP and Serpent2; as it will be emphasis in the following, the differences in these 2 examples is generally limited to 3 or 4 lines to change, whatever the complexity of the system. A README file suggests an order to execute examples. For each example, a compilation line is put at the end of the file ; this line is based on the following environment variable definitions as defined in section 1.1.1 ; theses variables should be set in your bashrc (or cshrc, ...).

Here, 2 examples are described ; other examples are described in this User Guide (see for example, § 4.8.1, 4.8.2, 4.8.3 and 4.5).

In MURE/mure3/share/examples directory, one finds static examples (i.e. without burn-up calculations) and their outputs7 (in MURE/examples/Static/Reference_Outputs/ directory). To run some of these examples you must provide an MCNP xsdir file, whereas for other ones user has to build the BaseSummary.dat and AvailableReactionChart.dat files. The “evolution” examples are in MURE/mure3/share/examples/Depletion/ directory. Again you will find an “Reference_Outputs” directory that contains directory output of the examples (also gzip). In these result directories, one has suppressed the MCNP “o” and “m” files but left the MCNP input files as well as the MURE results files.

1.1.6.1 basic MURE possibilities

This example (MURE/mure3/share/examples/Static/Putin_example.cxx and MURE/mure3/share/examples/Static/Putin_example_serpent.cxx) shows basic geometrical methods.

1.1.6.2 Fuel Evolution in a Sphere

The aim of this example8 of evolution (MURE/examples/Depletion/EvolvingSphere.cxx and MURE/examples/Depletion/EvolvingSphere_serpent.cxx) is to show basic burn-up calculation in MURE. The example describes a 50cm radius sphere with a inner sphere (R=30cm) of metal uranium, a middle ring (thickness=10cm) of UOx fuel and an external water reflector (thickness=10cm). The water will not evolve with time ; the thermal power of this “reactor” is 100 MW and it is kept constant for the whole evolution performed from t = 0 to t = 3 years.

1.2 MURE Package structure

The distribution package contains:

The general structure of main MURE’s classes is shown in figure 1.2.

pict

Figure 1.2: Main class structure of MURE. Main attributes are shown in blue. A Shape_ptr is only a typedef of Reference_ptr for a Shape. Idem for Nucleus_ptr.

As it can be seen in figure 1.2, the link between classes is assured by the main class MURE. The Shape_ptr and Nucleus_ptr of the figure are smart pointer on Shape and Nucleus objects. The figure 1.3 shows the Shape inheritance graph using the 2 name spaces associated to the MC transport code used (MCNP or Serpent).

pict

Figure 1.3: Shape inheritence scheme using with the 2 namespaces.

These name spaces are very important (see §3.1) ; one must use 1 (and only one) of them in any MURE input file. Results from MC run (Serpent detectors or MCNP tallies) are stored in the MureTally class which contains MureBin (cell, surface, ... bins) and TallyMultiplicator for reaction rates (see fig. 1.4 and 1.5).

pict

Figure 1.4: MureTally class.

pict

Figure 1.5: MureTallyBin class.

1.3 MURE class

The MURE class is some kind of super class that handles connections between all the other classes. It controls all flows: MC input file writing, nuclei trees building, evolution, ...

A global pointer (gMURE) on this class is defined to allow interaction between classes. For example, when the geometry has been defined, as well as the MC source, tallies (or detectors), ...., the MC file will be written using

gMURE->BuildMCFile("myfile");

this will generate an MC input file called “myfile” (if no argument is given to the previous method, the default input file name is “inp”).

The name of MC exec command is given by

gMURE->SetMCExec("mcnp");
or

gMURE->SetMCExec("sss2");

By default, MCNP::Connector set the the default value for the MC executable to “mcnp6” and Serpent::Connector set it to “sss2”.

MCNP Only: You can specify the particle running mode by either gMURE->SetModeN() (the default) to run MCNP with neutron transport, gMURE->SetModeNP() to transport neutrons and photons and so on. This is not take into account in the Serpent version.

Many methods exist and are grouped in sections to help the users ; therefore examination of MURE class header is strongly recommended (see Doxygen class description).

1.4 MURE basic files

Different files are used or built in MURE. It is important to know what these files are to understand the behavior of MURE.

1.4.1 Files of “MURE/mure3/doc”

This directory contains the most updated user guide and class documentations.

1.4.2 Files in “MURE/data”

In the root MURE tree, there is a very important directory called “data” ; this directory is copied by the buid.sh script to MURE/mure3/share/data (known in MURE as the DATADIR) ; this directory can be set or changed by using either the environment variable DATADIR

export DATADIR=”$MURE_DIR/share/data”

or using the MURE::SetDATADIR() method10. In this directory, you will find:

NOTE: When you add new nucleus or temperature in nuclear cross-sections files, these modifications will be taken into account ONLY if BaseSummary.dat and AvailableReactionChart.dat have been updated. You also need to remove the local ReactionList directory (see after § 1.4.7).

1.4.3 Files of “MURE/source/utils”

This directory contains codes that can be very helpful to users.

1.4.4 Graphical User Interface “MURE/source/gui”

The GUI of MURE is located in the MURE/source/gui directory. The executable name is “MureGui” in the MURE/mure3/bin directory (if you ask for it in the installation procedure) and it reads the result of a MURE or Dragon evolution. To use it, user must install ROOT (free download at https://root.cern.ch). If the Lapack package is needed, (see http://www.netlib.org/lapack/index.htmll)12. Using MureGui without argument gives a short description of the code used (see also § 8.3 ).

1.4.5 “MURE/source/external” and “MURE/source/thirdparty”

1.4.6 MURE Source files

The MURE source files are located in “MURE/source/include” for the header files (.hxx) and in “MURE/source/src” for implementation (.cxx). When user is modifying the MURE source, he has to recompile the MURE package (by a “build.sh” script or cmake13). In MURE/source/external, one finds 2 auxiliary libraries, ValErr and Mctal that are used in MURE. The ValErr defines a class to handle numbers and their errors, and Mctal is dedicated to reading/writing MCNP “m” files.

In any program using MURE, you have to include at least:

#include >MureHeaders.hxx>

and either one of the following

#include <MureMCNPHeaders.hxx>

or

#include <MureSerpentHeaders.hxx>

If needed on can add also :

#include <libmctal/TMTally.hxx> 
#include <libmctal/TMComment.hxx> 
#include <libmctal/TMctal.hxx> 

1.4.7 Other files

Each time a new evolution is run, a directory ReactionList is locally built (if it is not already there) with the available reaction of user’s nuclear data base (one binary file for each Z). A list of suppressed reactions for some of the nuclei (because they are lower than a given threshold) is also written in ReactionList/SuppressReaction.dat. The aim of theses files is to save time when using the same nuclear data base AND the same reaction threshold, life-time cutting and so on. Thus if you modify one of these, DO NOT FORGET to remove the ReactionList directory.

Chapter 2
What’s new in MURE

  1. January 2024:

    1. add a MureConsole, a kind of “text” version of MureGui to extract as text data from a single evolution directory ; The Console class, called by MureConsole, should be modified to extract your particular data ; it provided an example which should by easy to modify for your own purpose. The compilation is down via make -f Makefile.console.

    2. Modify the GenerateFPYield (version 2) to be able to generate FP independent yields from different ENDF libraries where data are provided as a single file or in multiple files ; all the files of the provided directory will be used as ENDF file format to extract FP yields for MURE evolution. A ReadFPyield code can be used to extract a given FP yield distribution as a function of the mass.

  2. December 2021: MURE/SMURE Version 3!

    1. Use smart pointer for gMURE and gTREE : each ZAI is now unique and only provided by the gTREE smart pointer.

    2. Evolution data can be written (when ASCII is chosen) with the JSON format via the EvolutionSolver::SetWriteJsonFormat().

    3. Complete MURE directory tree reorganization:

      • the root MURE tree contains: cmake, data, doc, examples, source ans tests directories as well as CMakeLists.txt and build.sh (for cmake install), LICENSE and README.md file.

      • source directory contains: src, include, utils, gui and external

      • the new install procedure is based on cmake (using either direct cmake command or the recommended build.sh install script)

      • and the default install directory is now mure3 (in MURE root directory) with a share dircetory (containing a copy of examples, data and doc (pdf, html)), include, lib and bin (with MureGui) ; this evolution allows a clear separation between MURE sources and user’s files.

      • change environment variables setting to take theses changes into account (see build.sh script for correct setting)

    4. Modernization of the code (all source code) according to c++11 standard.

    5. Many small (or not) bugs have been corrected (in particular in ReactorAssembly) and a great effort has been done in fixing memory leaks (rewrite of many constructors and destructors, used of smart pointers)

    6. documentation update.

  3. February 2021

    1. rewrite examples tree structure

    2. update doc (lyx & html & User Guide & doxygen)

    3. change compilation line in examples: it is now based on 3 environment variables (that should be defined in your bashrc or cshrc or ...) :

    export MURE_DIR="/path_2_MURE_install_dir/"
    export MURECXXFLAG="-g -fPIC -std=c++14 -I$MURE_DIR/source/include -I$MURE_DIR/source/external \
                        -I$MURE_DIR/source/thirdparty/eigen"
    export MURELIBFLAG="-L$MURE_DIR/lib -lMURE -lvalerr -lmctal"

  4. January 2021

    1. Implementation of the CRAM (Chebyshev Rational Approximation Method) method to solve Bateman equations (which is the new default to solve Bateman equation). This new implementation has been done by Maarten Becker.
      This implies some changes in method names and reorganization of Evolution related classes :

      1. RK discretization is now called sub-steps discretization and EvolutionControl methods using “RK” use now “SubStep”.

      2. Two new classes, CRAMSolver and RKSolver have been created and the depletion method is chosen by either

         gMURE->GetEvolutionSolver()->SetDepletionSolver("RungeKutta")  
      or  
         gMURE->GetEvolutionSolver()->SetDepletionSolver("CRAM").

      The latter is the default. These classes inherit from BatemanSolver and the previous DynamicalSystem class is removed.

    2. A new Timer class has been added to store time spent (to replace OMP omp_get_wtime(), ...).

    3. The compilation of SMURE require at least “c++11” standard.

    4. MureGui update: try to read DATA_* files if BDATA_* files are not found (avoid the use of “-type A”).

    5. The MURE library has been rename libMUREpkg.so as libMURE.so : you must change the linked library name in your project

    6. a cmake installation is now available (see paragraph 1.1.1)

    7. remove 32 bit systems support

  5. February 2019

    1. Add a new button in MureGui (Inventory Tab) to plot inventories in atom percent of the total of selected nuclei ; this allows for example to obtain the plutonium vector at any time of an evolution.

  6. November 2016:

    1. Add Z revolution axis Torus for MCNP and Serpent ; see MathZTorus , MCNPZTorus and SerpentZTorus classes.

    2. Update “chart.jef3T” file that contains periods and decay modes for all the chart ; this file (JEFF 2 lib.) is kept unchanged, but a new file has been created “chart.jeff3.1.1” (from ref. [19]); it is the new default.

  7. November 2015:

    1. ReactorMesh and ReactorChannel class have been removed ; a new GenericReactorAssembly, and its concrete versions MCNP::ReactorAssembly and Serpent::ReactorAssembly replace the ReactorMesh class. It is a complete rewriting of that class, bu the spirit is kept (see § 9.1 and 9.3.2). It adds lot of new possibilities (add duct to an assembly, filling core with a ReactorAssembly, ...).

    2. COBRA class is renamed COBRA_EN and rewritten ; again backward compatibility is broken (it does not inherit from GenericReactorAssembly) but the spirit is not changed.

  8. October 2015:

    1. Correct a bug when reading the “standard” xsdir of MCNP5 ; isomeric cross-section libraries are now in accordance with the MCNP5 documentation.

  9. July 2015:

    • New MURE version 2.0 (SMURE). This is a major release of MURE. It couples now both MCNP & Serpent. Due to technical, philosophical difficulty and also to beauty in programming (dry codes), backward compatibility is broken. But the changes to user input files between MURE v1 and MURE v2 are small. See 3.

  10. November 2014:

    1. Add the possibility to use standard MCNP tallies 238U and multi-group tallies for other nuclei (see 5).

  11. October 2014:

    1. Add an optional equilibrium treatment for 135Xe (see 7.7)

  12. April 2013:

    1. Version of base and isomeric states (metastable and ground state) for 242Am in BaseSummary.dat file for standard library (distributed with MCNP or available@NEA) are now more conform to the reality ; you have to regenerate your BaseSummary.dat/AvailableReaction.dat files with the exec of MURE/utils/datadir. DON’T FORGET TO REMOVE YOUR ReactionList DIRECTORY IN ORDER TO TAKE INTO ACCOUNT THESE MODIFICATION.

    2. BasePriority has been debug and it seems to work as expected...

  13. March 2013:

    1. Add a generalization of the MURE::SetMode() method to allow any type of particle transport (mainly for MCNPX)

    2. Cell importance may have different value for each transported particles (by successive call to Cell::AddParticle method)

    3. extend the MCNPSource possibilities (particle distributions are now allowed for MCNPX)

  14. July 2012:

    1. Generalization of isomer production by (n,γ) for special cases (such as 241Am, 109Ag, …), and bug correction in (242Cm, 242Pu) production, see 6.1.6.

    2. Improve MurGui Spectrum radiotoxicity tab

    3. Allow to plot neutron balance in the “Reaction Rate” tab.

  15. June 2012:

    1. Add functions in MureGui (see MureGui’s Radiotoxicity tab)

      1. Nuclei extraction is now possible after a given cooling time.

      2. γ,β, α and neutron spectra for evolving materials can be computed.

    2. New MURE install procedure (more simple, more robust): any fresh install or update NEED to run first the install.sh script (require bash shell).

    3. Add the possibility to used multi-threading (OpenMP) during the evolution and for updating multigroup σϕ when the gcc version support this option (see INSTALL procedure and install.sh script and section 7.2.2.2).

    4. Add some new kind of MCNPSource (see section 5.2.3)

    5. Add fluence to dose conversion for tallies

  16. May 2011:

    1. Add a new class GammaSpectrum class (section 8.4.2)

    2. New option of MureGui : use the “radiotoxicity tab” of MureGui on a dumped Material created with MURE (see 8.3.6 and “MURE/example/GammaSpectrumExample.cxx”)

    3. New methods MCNPSource::SetAXS and MCNPSource::SetVec allows user to define collimated sources.

    4. In Tally::Tally(int type,const char *particle) : if type < 0 , it’s change the Tally units.

      1. For example Tally *t=new Tally(-4,"P") define a tally in MeV∕cm2 rather than Particles∕cm2 .

    5. A new method MURE::SetModeP() can be used to allow only the photon transport.

  17. July 2009:

    1. Improve the English in User Guide: many thanks to Erica Agostinho for this painful work.

    2. Implementation of PseudoMaterial: in order to take into account temperature effects, one can process the nuclear data at the desired temperature (using NJOY ) or use an interpolation between 2 already existing temperatures ; this later method is used in the PseudoMaterial techniques.

  18. Avril 2009: MURE is available at NEA DataBank

  19. January 2009:

    1. Improve reading time of BaseSummary.dat file : it is now greatly recommended that this file is ordered by Z,A,I (this is the case if it is generated by ExtractedXsdir.cxx and ExtractTree.cxx).

    2. Rewrite long parts of the documentation (User Guide and HTML class description)

    3. Modify examples directory: it now contains documented examples and output.

  20. November 2008: Improve radiotoxicity post treatment in MureGui.

  21. September 2008: Implement “multigroup” calculation for reaction rates: a very narrow energy binning for flux calculation is put in each MCNP run ; then reaction rates are obtained by reading ACE MCNP files after each MCNP run. The result of such method saves a considerable amount of CPU time for MCNP (at least a factor 30) with only a low percentage of discrepancy (∼1 to 3%) in result compared to the standard calculation (reaction rates are tallied in MCNP).

  22. April 2008:

    1. implementation of Predictor-Corrector method in the evolution

    2. possibility to read ASCII nuclear data file (ACE format): this avoids problems due to binary compatibility (size of real (float or double?), little or big endian, ...). BUT it is much longer to (1) build the ReactionList directory and (2) run a MCNP.

  23. December 2007:

    1. Disable the σϕ extrapolation by default. It has been shown, but not really understood, that this treatment introduces a larger dispersion in the result, after N identical evolutions, than doing nothing (no σϕ extrapolation).

    2. Modify the evolution using a MCNP User Input geometry file:

      1. “like but” cells can be used (but not evolved)

      2. “MCNP Transformation” cards can be used

      3. the 3rd block of the MCNP file is read and copied except the materials (that must be defined in the MURE file). Thus the source as well as all other cards of this block can be used without defining them in the MURE file (this is also true for user defined tallies).

  24. October 2007:

    1. Switch from ccdoc to Doxygen for class documentation

    2. Change completely the MURE directory trees

    3. NON BACKWARD COMPATIBILITY

      1. Material definition has changed: now only 2 constructors must be used:

        • Material(): for standard material (you have to specify density, ... with the Set methods)

        • Material(int): for using materials from an MCNP input file geometry.

      2. THE COPY CONSTRUCTOR HAS NOT TO BE CALL: CALL only the Material::Copy().

      3. The Material::Mix has been modified (number of arguments, units required). see Material.hxx

      4. The Material::AddNucleus has been modified (number of arguments). By default the proportion units are "kpMOL" (i.e. molar proportion). But the unit must be specified if you use a moderator (MT card of MCNP).

      5. The proportion units (both for Density and Proportion) must be used for any Material::GetProportion() and Material::GetDensity() methods

        • for Proportion the only valid units are: kpMOL(molar prop), kpMASS (mass prop), kpATCM3 (at/cm3), kpATBCM(at/barn.cm)

        • for density the only valid units are: kdGCM3 (g/cm3), kdATCM3 (at/cm3), kdATBCM (at/barn.cm)

      6. A new material has been defined : ControlMaterial (public of Material). This class is used for Poison, Fissile or other control of reactivity (e.g. poison.cxx in MURE/examples)

      7. All Print() methods now return a string instead of a void: to use them: cout<<Mat->Print()<<endl; for example.

      8. EvolutionControl class has been cleaned (as well as MURE class). If you want to use special control you have to write your own derivative class. Examples using PoisonControl, FissileControl & HNControl (but Adrien you have to rewrite them and look carefully at TMSR.cxx in MURE/examples.) and Rod control are defined in source/src. You can used them as they are or defined your own using these examples.

    4. Almost all cout/cerr have been removed from classes; used instead LOG_DEBUG, LOG_INFO, ... ; LOG_INFO is now independent of LogLevelMessage ; it is always printed. If MURE::SetMessageLevel is set to LOG_LEVEL_DEBUG, all LOG_DEBUG are printed. But if MURE::SetSilentDebug is used, only the LOG_DEBUGs of methods where a "int DODEBUG=1" is inserted are printed. Thus, using LOG_DEBUG, avoids to comment all "cout" when no debugging is desired.

    5. ALL EXAMPLES HAVE BEEN UPDATED TO TAKE INTO ACCOUNT THESE MODIFICATIONS: please READ THEM!!!!!!!!!!

  25. September 2007:

    1. Correct an important bug in evolution using an MCNP User Input geometry file: number of Materials were not correct (Thanks to Jan Frybort).

  26. Juin 2007:

    1. Rename the MURE header file Shapes.hxx into MureHeaders.hxx : this is more logical...but you have to change your MURE files....

    2. Add a new class EvolutionWrapper to simplify and extend EvolutionControl capabilities

    3. Suppress the writing of BDATA_xxx and DATA_xxx ; now, by default, only BDATA_xxx are written. This can be changed using the MURE::SetWriteBinaryData() and MURE::SetWriteASCIIData() methods.

    4. One can start an evolution from a given step : suppose that the evolution stops at step i ; an evolution can be started from the step i+1 using MURE::Evolution(T,i+1). Warning: it is probably not correct for OutCoreEvolutiveSystemVector...you must do the evolution from the first step as before.

Chapter 3
From MURE to SMURE : Choosing the Monte-Carlo Transport Code

As already said, coupling MURE with a other MC code had 2 consequences:

The aim is with a minor change to input file, it will be able to generate output for MCNP or for Serpent 2. In this chapter, it is explained first how to switch from an MCNP input file to a Serpent input file (or vice-versa) and then how to migrate from MURE 1.0 to SMURE (MURE 2.0) .

3.1 Switching from MCNP to Serpent in a MURE input file

SMURE has been designed to make the change from Serpent input to MCNP input (or vice-versa) as simple as possible. Of course, Serpent is not MCNP and conversely. Thus each of these MC codes have their own possibility, not implemented in the other. Moreover, the maturity of the MURE coupling is much greater (few years of work for MCNP versus 4 months of work for Serpent). Thus this section deals only with things that are possible with both MC codes and has been implemented in MURE. Nevertheless, it seems to us that is it enough to make a standard reactor study (including evolution).

3.1.1 Principle of the implementation

In MURE, the choice of MC output is performed via a ConnectorPlugin ; this class has 2 sibling, a MCNP::Connector and a Serpent::Connector (see Fig. 3.1) ;

pict

Figure 3.1: The ConnectorPlugin class.

thus the main trick is the use of a specific namespace (MCNP or Serpent). Geometrical shapes are defined in MURE as MathShape (i.e. a “pure” mathematical Shape, e.g. a MathSphere). Then, in each namespace, are defined concrete (or real) shapes that essentially provide a print method derivating from the MathShapes (see Fig. 1.3). For example, a Serpent::Sphere and MCNP::Sphere are daughters of MathSphere and implement a Serpent::Sphere::Print() or a MCNP::Sphere::Print() that will produce respectively the right output for the corresponding MC code. The same trick applies for tallies (or detectors in Serpent).

Thus, using a Connector, a Sphere, ... in MURE is just “resolved” by the namespace (that should not be forgotten). In order it works without problems, appropriate header files must be included. First header file, is the “pure” MURE header : MureHeaders.hxx which includes all pure MURE classes. It is thus included in all MURE input files, whatever will be the output. The second header file is linked to the output ; MureMCNPHeaders.hxx and MureSerpentHeaders.hxx includes all classes that couple MURE with MCNP and Serpent respectively. They include particularly MCNP::Connector (Serpent::Connector), MCNP::Sphere (Serpent::Sphere), ...

Of course connectors, because of some profound differences between MCNP & Serpent, have specific methods existing only in one type of connector. The same applies also for the Tally class which is more complete in MCNP. The last case is for particle sources which are very different from MCNP to Serpent. This is why, sibling of MCSource class keep the MC code prefix in there name (MCNPSource and SerpentSource) in order that the user don’t forget its specificity. An other cause of problems, switching from one code to the other is the MURE::AddSpecialCard() which allows to use specific Serpent or MCNP cards that have not been implemented in MURE. Thus a careful check of Connector, MCSource and MURE::AddSpecialCard() should be done when switching from one code to the other.

3.1.2 A example

Switching from one MC code to the other can be performed very fast : edit the file (let say a MCNP version) and replace all “MCNP” string by “Serpent” string. This should change 4 occurrences: the header, the namespace and the source name (See example of Tab. 3.1). Then remove specific command for MCNP and add specific command for Serpent. To illustrate this, let us consider the “evolving sphere” example. The “red” parts of Table 3.1 correspond to a “replace all” “MCNP” string by “Serpent” string. The magenta words are “output” directory name, that here, are different to keep a trace of both evolution. The blue lines is an example of a pure specific method of one code (here MCNP) of one of the 3 classes Connector, (MCNP/Serpent)Source or Tally. And then, in green, a MURE::AddSpecialCard() that in this case is only for Serpent.

Table 3.1: Switching from MCNP (left) to Serpent (right) output files.

#include <iostream>
using namespace std;
#include <MureHeaders.hxx>
#include <MureMCNPHeaders.hxx>
using namespace MCNP;
int main(int argc, char **argv)
{
  Connector *plugin=new Connector();
  //remove the MCNP "r" file after each evo step
  plugin->SetRemove_r_files();
  gMURE->SetConnectorPlugin(plugin);
//
//  Nuclei tree construction & nuclear data library
//
  gMURE->SetMessageLevel(3);
  gMURE->SetDATADIR("../../data/");
  gMURE->SetAutoXSDIR();
  gMURE->SetShortestHalfLife(3600);
  gMURE->GetTemperatureMap()->SetDeltaTPrecision(1500);
  gMURE->KeepOnlyFissionProductSelection();
//
//  Materials
//
  Material *H2O= new Material();
  H2O->SetDensity(1.);// density in g/cm3
  H2O->AddNucleus(1, 1, 0, 2.); // H of H2O
  H2O->AddNucleus(8, 16, 0, 1.);// O of H2O
  
  Material *Uranium_metal= new Material(); 
  Uranium_metal->SetDensity(19.);          
  Uranium_metal->SetTemperature(600.); //Temperature in K
  Uranium_metal->AddNucleus(92,235,0.20);//U-235 in mol
  Uranium_metal->AddNucleus(92,238,0.80);//U-238 in mol
  Uranium_metal->SetEvolution(); //this is an evolving material    
  
    Material *UOx= new Material(); 
  UOx->SetDensity(10.2);
  UOx->AddNucleus(92,235,0, 0.05);// 5%(in mol) of U-235
  UOx->AddNucleus(92,238,0, 0.95);// 95%(in mol) of U-238
  UOx->AddNucleus(8,16,0, 2.);  // 2 O for 1 U02
  UOx->SetEvolution();  //this is an evolving material    
//
//  Shapes
//
  Shape_ptr Shape_MySphere(new Sphere(0.5));
  Shape_ptr Shape_Exterior(!Shape_MySphere);
  Shape_ptr Shape_MySphere2(new Sphere(0.4));
  Shape_ptr Shape_MySphere3(new Sphere(0.3));
 
  Shape_MySphere3>>Shape_MySphere2>>Shape_MySphere;

//
//  Cells
//
  Cell* CellSphere  = new Cell(Shape_MySphere , H2O);
                            
  Cell* CellSphere2 = new Cell(Shape_MySphere2, UOx);
 
  Cell* CellSphere3 = new Cell(Shape_MySphere3, Uranium_metal);
 
  Cell* CellExterior = new Cell(Shape_Exterior, 0, 0);
//
//  General Cards
//
  //define the neutron source (KCODE)
  MCNPSource *s=new MCNPSource(500);
  s->SetKcode(60,30,1.);
    s->SetEnergy(2e6);
  gMURE->SetSource(s);
  
  //  MC option
  gMURE->SetComment("Evolution of a Sphere");
  gMURE->SetMCInputFileName("inp");
  gMURE->SetMCRunDirectory("sphere","keep");
  
//
// Evolution
//
  gMURE->SetPower(1e8);  //give the total power: 100 MW
  
  double TEnd=3*365.25*3600.*24;
  
  vector<double> T;
  int Nt=6;
  double dt=TEnd/Nt;
  
  for(int i=0; i<Nt;i++)
      T.push_back(i*dt);
 
  gMURE->Evolution(T);
}

#include <iostream>
using namespace std;
#include <MureHeaders.hxx>
#include <MureSerpentHeaders.hxx>
using namespace Serpent;
int main(int argc, char **argv)
{
  Connector *plugin=new Connector();
  
  
  gMURE->SetConnectorPlugin(plugin);
//
//  Nuclei tree construction & nuclear data library
//
  gMURE->SetMessageLevel(3);
  gMURE->SetDATADIR("../../data/");
  gMURE->SetAutoXSDIR();
  gMURE->SetShortestHalfLife(3600);
  gMURE->GetTemperatureMap()->SetDeltaTPrecision(1500);
  gMURE->KeepOnlyFissionProductSelection();
//
//  Materials
//
  Material *H2O= new Material();
  H2O->SetDensity(1.);// density in g/cm3
  H2O->AddNucleus(1, 1, 0, 2.); // H of H2O
  H2O->AddNucleus(8, 16, 0, 1.);// O of H2O
  
  Material *Uranium_metal= new Material(); 
  Uranium_metal->SetDensity(19.);          
  Uranium_metal->SetTemperature(600.); //Temperature in K
  Uranium_metal->AddNucleus(92,235,0.20);//U-235 in mol
  Uranium_metal->AddNucleus(92,238,0.80);//U-238 in mol
  Uranium_metal->SetEvolution(); //this is an evolving material    
  
    Material *UOx= new Material(); 
  UOx->SetDensity(10.2);
  UOx->AddNucleus(92,235,0, 0.05);// 5%(in mol) of U-235
  UOx->AddNucleus(92,238,0, 0.95);// 95%(in mol) of U-238
  UOx->AddNucleus(8,16,0, 2.);  // 2 O for 1 U02
  UOx->SetEvolution();  //this is an evolving material    
//
//  Shapes
//
  Shape_ptr Shape_MySphere(new Sphere(0.5));
  Shape_ptr Shape_Exterior(!Shape_MySphere);
  Shape_ptr Shape_MySphere2(new Sphere(0.4));
  Shape_ptr Shape_MySphere3(new Sphere(0.3));
 
  Shape_MySphere3>>Shape_MySphere2>>Shape_MySphere;

//
//  Cells
//
  Cell* CellSphere  = new Cell(Shape_MySphere , H2O);
                            
  Cell* CellSphere2 = new Cell(Shape_MySphere2, UOx);
 
  Cell* CellSphere3 = new Cell(Shape_MySphere3, Uranium_metal);
 
  Cell* CellExterior = new Cell(Shape_Exterior, 0, 0);
//
//  General Cards
//
  //define the neutron source (KCODE)
  SerpentSource *s=new SerpentSource(500);
  s->SetKcode(60,30,1.);
  s->SetEnergy(2e6);
  gMURE->SetSource(s);
  
  //  MC option
  gMURE->SetComment("Evolution of a Sphere");
  gMURE->SetMCInputFileName("inp");
  gMURE->SetMCRunDirectory("sphere_serpent","keep");
  gMURE->AddSpecialCard("plot 3 2000 2000");
  
//
// Evolution
//
  gMURE->SetPower(1e8);  //give the total power: 100 MW
  
  double TEnd=3*365.25*3600.*24;
  
  vector<double> T;
  int Nt=6;
  double dt=TEnd/Nt;
  
  for(int i=0; i<Nt;i++)
      T.push_back(i*dt);
 
  gMURE->Evolution(T);
}

As Shown on this example, the number of line to modify is very limited ; moreover, dealing withe a more complex case, this number will not increase significantly (probably less than a factor 2).

3.2 How migrate my old MURE V1.x file to the MURE V2.x file

3.2.1 What has definitely change

In this part, only the coupling with MCNP is considered (as in the old MURE v1). Main changes are due to

Table 3.2: Summary of changes from MURE to SMURE.

MURE version 1

MURE version 2 : SMURE

#include <MureHeaders.hxx>

#include "MureHeaders.hxx"

#include "MureMCNPHeaders.hxx"

using namespace MCNP;

int main(...)

int main(...)

{

{

  Connector *plugin=new Connector();

  gMURE->SetConnectorPlugin(plugin);

  gMURE->SetRemove_r_files();

  plugin->SetRemove_r_files();

  gMURE->BuildMCNPFile();

  gMURE->BuildMCFile();

  gMURE->SetMCNPInputFileName("inp");

  MURE->SetMCInputFileName("inp");

  gMURE->SetMCNPRunDirectory("sphere","keep");

  gMURE->SetMCRunDirectory("sphere","keep");

  gMURE->SetMCNPExec("mcnp5");

  gMURE->SetMCExec("mcnp5");

  gMURE->SetWriteASCIIData();

  gMURE->GetEvolutionSolver()->SetWriteASCIIData();

  gMURE->UseMultiGroupTallies();

  gMURE->GetEvolutionSolver()->UseMultiGroupTallies();

  gMURE->SetUseNewDBCN();

  gMURE->SetUseNewRandomSeed();

  Cell *MyLattice=new Cell(AGenerator);

  LatticeCell *MyLattice=new LatticeCell(AGenerator);

  MyLattice->Lattice(1,-rx,rx,-ry,ry,-rz,rz);

  MyLattice->SetLatticeRange(-rx,rx,-ry,ry,-rz,rz);

}

}

The time spent to change a very big file from version 1 to version 2 is around 5 minutes maximum.

3.2.2 What is still existing but can be done in a more elegant way

In most of real case fro reactor physics, one has to use lattice, fuel pin, ... A new improvement of MURE version 2 is to implement a PinCell which in fact copy the elegant Serpent’s pin definition. To illustrate the advantage of such a PinCell class, the main part of the lattice example in the “old” way is presented and then rewrite using PinCell. For sake of clarity, only Shape and Cell part are presented. All sizes have been defined, as well as 3 Materials : H2O, Iron, and UOx. As shown on Tab. 3.3, the use of PinCell makes the code shorter (the 8 blue lines of the shape declaration have been suppressed), and easier to define and to read. The PinCell declaration make cell part slightly longer (4 ”magenta” lines versus 6 “green” lines). To be noticed, the use of a PinCell without a layer allows to define the whole universe. And a last comment on the FillLattice methods : in the first version, a Shape is used as argument (red variables on the left) where as in the new version this is a PinCell (red variables on the right).

Table 3.3: Comparison of a lattice in the Old Mure version1 (left) and with the PinCell of SMURE(right).

// 
//Shapes 
// 
//the Vessel : a full Tube 
Shape_ptr Vessel(new Tube(VesselH/2, VesselR)); 
//the Exterior 
Shape_ptr Exterior(!Vessel); 
//the lattice generator 
Shape_ptr LatticeGenerator(new Hexagon(HexaH/2,HexaSide)); 
LatticeGenerator->SetUniverse(); 
 
// put the lattice generator in the vessel 
LatticeGenerator>>Vessel; 
 
// a small cylinder that will contain UOx 
Shape_ptr Pin(new Cylinder(PinR)); 
Pin->SetUniverse(); 
// a bigger one that define fuel clad 
Shape_ptr Clad(new Cylinder(CladR)); 
Clad->SetUniverse(Pin->GetUniverse()); 
//Outside of the CladShape_ptr OutClad=!Clad; 
Pin>>Clad; 
//the whole space that is use near the vessel border 
Shape_ptr Whole(new Node(-1)); 
Whole->SetUniverse(); 
// 
//Cells 
// 
 
Cell *pin=new Cell(Pin,UOx); 
Cell *clad=new Cell(Clad,Iron); 
Cell *surrounding_clad=new Cell(OutClad,H2O); 
 
Cell *whole=new Cell(Whole,H2O); 
 
//Exterior of the vessel: void and importance=0 
Cell *exterior=new Cell(Exterior,0,0); 
Cell *vessel=new Cell(Vessel); 
vessel->SetComment("The vessel"); 
//the lattice 
LatticeCell *Pavage=new LatticeCell(LatticeGenerator); 
Pavage->SetLatticeRange(-range,range,-range,range); 
for(int i=-range; i<=range; i++) 
for(int j=-range; j<=range; j++) 
{ 
int pos[3]={i,j,0}; 
double xt = HexaSide*sqrt(3.) * (i + j*0.5 ) ; 
double yt = 1.5 * HexaSide * j; 
double X[2]={xt,yt}; 
if(IsHexagonInTube(X,LatticeGenerator,0.9)) 
Pavage->FillLattice(Pin,pos); 
else 
Pavage->FillLattice(Whole,pos); 
}

// 
//Shapes 
// 
//the Vessel : a full Tube 
Shape_ptr Vessel(new Tube(VesselH/2, VesselR)); 
//the Exterior 
Shape_ptr Exterior(!Vessel); 
//the lattice generator 
Shape_ptr LatticeGenerator(new Hexagon(HexaH/2,HexaSide)); 
LatticeGenerator->SetUniverse(); 
 
// put the lattice generator in the vessel 
LatticeGenerator>>Vessel; 
 
 
 
 
 
 
 
 
 
 
 
 
// 
//Cells 
// 
PinCell *pin=new PinCell; 
pin->AddLayer(UOx,PinR); //layer 0 = fuel 
pin->AddLayer(Iron,CladR); //layer 1 = the clad 
pin->SetSurroundingMaterial(H2O);//surrounding the clad = H2O 
PinCell *whole=new PinCell; 
whole->SetSurroundingMaterial(H2O); //no layer in that PinCell 
 
//Exterior of the vessel: void and importance=0 
Cell *exterior=new Cell(Exterior,0,0); 
Cell *vessel=new Cell(Vessel); 
vessel->SetComment("The vessel"); 
//the lattice 
LatticeCell *Pavage=new LatticeCell(LatticeGenerator); 
Pavage->SetLatticeRange(-range,range,-range,range); 
for(int i=-range; i<=range; i++) 
for(int j=-range; j<=range; j++) 
{ 
int pos[3]={i,j,0}; 
double xt = HexaSide*sqrt(3.) * (i + j*0.5 ) ; 
double yt = 1.5 * HexaSide * j; 
double X[2]={xt,yt}; 
if(IsHexagonInTube(X,LatticeGenerator,0.9)) 
Pavage->FillLattice(pin,pos); 
else 
Pavage->FillLattice(whole,pos); 
}

Thus user is strongly encourage to use these PinCells.

Chapter 4
Geometry Definition

4.1 Introduction

C++ logic for building geometries is slightly different for the MCNP or Serpent one ; therefore, each time a new geometry is built you should check it with MCNP/Serpent before using it.

There are two base classes to build a geometry: Shape and Cell. Shape describes only geometrical shapes, and Cell corresponds to an MC code cell (i.e., it has a material, importance, etc.).

Shape objects correspond to simple geometrical shapes (sphere, plane, ...) as well as more complex ones resulting from the intersection and/or union of simple shapes (Intersections/unions are defined by the Node class). A Node is a “tree” of intersections/unions of Shapes. For fast calculations, a node tree has to be as simple as possible. Special methods are available for simplifying the node trees which can (in general) determine whether or not a Shape is disjointed (or included) of (in) another Shape (see example in figures B.1 to B.4 in Appendix A.6).

These simplifications may result in the deletion of some Shapes. But, because one must not destroy a Shape in a Node if that Shape belongs to another Node, a special way to handle Shape creation/destruction has been implemented (via Reference_ptr, which is a C++ smart pointer).

In conclusion : the user must only use Shape_ptr (a “reference Shape”) and not Shape. Shape_ptr is a pointer to Shape with Reference_ptr. In that way, the Shape_ptr will be destructed only if it is no longer referenced, otherwise, its “deletion” leads to a decrement of the number of references.

4.2 Definition of geometrical shapes

A Brick is a rectangular parallelepiped (and can be infinite). A Tube is a finite cylinder with an optional inner radius that defines a “tube”. Hexagons can have finite or infinite height. In general, a user will only need to use Spheres, Bricks, Tubes and Hexagons.

One can define the interior or the exterior of each Shape. The complement of a Shape can be defined by the Not() method as well as by the “!” operator (see examples given later on).

Other shapes available in MCNP(X) or Serpent may be added by any user. This is not too hard even if it requires some work. The best way to do it, is to read the implementation of existing shapes to have an idea.

4.2.1 Units

WARNING: In MURE, the length unit is the meter (whereas in MCNP and Serpent it is the cm). In particular, volumes in MURE are in m3.
Length m
Energy eV
Temperature K
Density g/cm3, atom/cm3 or atom/(barn.cm)
Proportion %mol, %mass, atom/cm3 or atom/(barn.cm)

4.2.2 Examples of simple shapes

To define the interior of a origin center sphere of radius R:

Shape_ptr S(new Sphere(R));

To define the outer part of S one can do

Shape_ptr Ext_S(!S);

or

Shape_ptr Ext_S(S->Not());

or

Shape_ptr Ext_S(new Sphere(R,0,0,0,1);

where 3 zeros correspond to the sphere center (the origin) and +1 to the exterior of the Sphere (default=-1).

For Hexagons and Bricks, two versions exist: finite or infinite Shapes. For finite bricks and hexagons, you should define it as:

Shape_ptr B(new Brick(HalfX, HalfY, HalfZ, Signe));
Shape_ptr H(new Hexagon(HalfHeight, Side, Signe));

where HalfX (resp. Y and Z) is the Half length of the brick in the X (resp. Y and Z) direction, HalfHeight the half height (!) of the hexagon, Side its side, and Signe the sign defining whether it is the inner or the outer part, just as for the sphere. For infinite ones, to avoid conflicts between the definitions, a string at the beginning is necessary (and, of course, HalfHeight and HalfZ are irrelevant) :

Shape_ptr B(new Brick(“any string you want”, HalfX, HalfY, Signe));
Shape_ptr H(new Hexagon(“any string you want”, Side, Signe));

The infinite shapes obtained are parallel to the Z axis, but they may be rotated afterwards.

4.2.3 Examples of simple Nodes

A Node is the Union (+1) or the Intersection (-1) of Shapes (i.e., simple Shapes or Nodes).

To define the intersection of a sphere centered at (x0,y0,z0) of radius R with a square brick of side a centered at the origin, one may use:

Shape_ptr S(new Sphere(R,x0,y0,z0);
Shape_ptr B(new Brick(a/2,a/2,a/2));
Shape_ptr Inter=S & B;

whereas the union of the sphere and the brick may be defined as

Shape_ptr Union=S | B;

4.2.4 Moving a Shape

It is possible to move (translation/rotation) a Shape (or a Node) via the Shape::Translate and Shape::Rotate methods. For example, to translate the Shape_ptr B of (dx,dy,dz),

B->Translate(dx,dy,dz);

and to rotate clockwise B around (x0,y0,z0) of ϕ, 𝜃 and ψ around the z, y and x axis respectively:

double center[3]={x0,y0,z0};
B->Rotate(phi,theta,psi,center);

Note that the translation/rotation of a Shape translates/rotates also the inside shapes (see next section).

To be noticed: Angles are in radians.

4.3 The “put in” operator ">>"

It is possible to put a Shape_ptr A inside another B one, via the “put in”1 operator A>>B ; this operator works differently depending on A:

Example:

Shape_ptr S(new Sphere(R));
Shape_ptr B(new Brick(a/2,a/2,a/2));
B>>S;

In this example, the Brick B is put inside S.

Note that now S and B are new Shape_ptr: S=S & !B (i.e. the sphere without the brick) and B=S & B (i.e, the intersection of the brick and the sphere) as already mentioned above (see fig. 4.1).

Note:

A>>B>>C; 

means that A is put in B and the result is also put in C. A is an inside shape of B and C ; B is an inside shape of C.

C>>D; 

Inside shapes of C are cleared and D has A, B and C as Inside Shapes.

4.4 Clone of a Shape

It may be sometimes very useful to clone a Shape, i.e., to create a new Shape with the same properties as the original one. Here is an example:

Shape_ptr B1(new Tube(1,0.5));
Shape_ptr T1(new Tube(1,0.4));
Shape_ptr T2(new Tube(1,0.3));
T2>>T1>>B1;
Shape_ptr B2=B1->Clone();
B2->Translate(1,0,0);

B2 has exactly the same aspect as B1 (with 2 tubes inside) but it is translated by 1m in the x direction, whereas B1 is centered at the origin.

4.5 Boundary conditions

In MCNP, it is possible to use special boundary conditions: mirror, white or periodic boundary conditions. In Serpent, boundary conditions can be applied only on the shapes associated to cells with the “outside” Serpent’s material. The only valid possibilities for Serpent are black (totally absorbing, the default), mirror or periodic.

To be noticed also, the “outside” cells of Serpent don’t support unions. For any geometry, if the outside cell (Serpent namespace only) is not a “simple” shape (simple in the Serpent’s sense, it includes Spheres, unrotated Tubes, Bricks or Hexagons), then the user define outside shape (the one which is in the 0 importance cell) is put into a Sphere.

In Serpent, boundary conditions only apply to Bricks (square, cuboid, ...) or Hexagons (hexxc, hexxprism, ...). Mirror and Black boundary conditions have been implemented in MURE for unrotated Bricks or Hexagons in Serpent. In any case, the use of boundary conditions in MURE could be used but have not been fully tested...so it is advised to use them carefully. In all cases, the exterior of the Shape_ptr with boundary conditions must have a zero importance. You can apply only one type of boundary conditions to a Shape_ptr. Here is a brief description of these boundary conditions:

4.6 Recommendations

4.7 Definition of MC cells

A Cell is defined by a Shape_ptr, and if needed a Material.

MC cells are defined via the Cell class ; one has to give

For example to construct a cell composed of a full tube of radius R and height H, made of B4C

Shape_ptr C(new Tube(H/2,R));
Cell* c=new Cell(C,B4C);

where B4C is a (Material*)3. In this example, one can define the exterior cell as

Shape_ptr Exterior(!C);
Cell* exterior=new Cell(Exterior,0,0);

The first zero means that the material is void and the second one means that neutron importance is set to 0 in the exterior cell (“outside” material for Serpent).

NOTE: If not only neutrons are transported in MCNP, it is useful to define the desired mode (i.e., NP, NE, NPE for Neutron and Photon, Neutron and Electron and Neutron, Photon and Electron) before any Cell definition ; indeed, in this case, the given cell importance applies to all particle types. For example:

gMURE->SetModeNP();
Cell *c=new Cell(C,B4C);
Cell *exterior=new Cell(Exterior,0,0);

will set the transport mode to neutrons and photons and the cell c will have an IMP:N=1 and IMP:P=1 whereas the cell exterior will have an IMP:N=0 and IMP:P=0.

4.7.1 Cell and clone shapes

Suppose we have defined 2 tubes B1 and B2 as follows:

Shape_ptr B1(new Tube(1,0.5));
Shape_ptr T1(new Tube(1,0.4));
Shape_ptr T2(new Tube(1,0.3));
T2>>T1>>B1;
Shape_ptr B2=B1->Clone();
B2->Translate(1,0,0);

the cell definition for B1 is not a problem:

Cell *b1=new Cell(B1,Graphite);
Cell *t1=new Cell(T1,Iron);
Cell *t2=new Cell(T2);

where Graphite and Iron are two Material*. Here, the result will be a tube of graphite containing a tube of iron containing a void tube. To define analog cells for B2 we have to do:

Cell *b2=new Cell(B2,Graphite);
Cell *tt1=new Cell(B2->GetOriginalInsideShape(1),Iron);
Cell *tt2=new Cell(B2->GetOriginalInsideShape(0));

Indeed, the “Inside shapes” of B2 are ordered from the most inner (B2->GetOriginalInsideShape(0) clone of T2), to the most outer (B2->GetOriginalInsideShape(1) clone of T1).

4.8 Lattice

Lattices are used in MC to fill cells of repeated structure. In MCNP, there are 2 types of lattice, hexahedra (type=1) or hexagonal (type=2). In Serpent, circular lattice are also defined but they are yet not implemented in MURE.

In this section, four examples of lattices are presented (from the most simple to the most complex case). Of course the lattice type does not change the declaration (except that the lattice generator is a Brick for hexahedra lattice whereas it is an Hexagon for the the hexagonal one). Lattice philosophy is somehow different from MCNP and Serpent ; the MURE implementation is closer to the MCNP’s one.

The general philosophy for lattice declaration is the following:

  1. In the Shape section

    1. define the Shape_ptr C that will be filled by the lattice (like core),

    2. define the Shape_ptr G that is used as lattice generator (a Brick or an Hexagon) and give to that generator a universe number (via Shape::SetUniverse()),

    3. put G in C with G>>C

    4. define all Shape_ptr Bi that shall be used in the lattice and assign them a universe number (via Shape::SetUniverse()) or use PinCell (see (a) of Cell Section)

  2. In the Cell section

    1. define PinCell Pi (and/or use (d) of previous part)

    2. build the Cell for Shape_ptr C

    3. build the LatticeCell with the generator Shape_ptr G

    4. [optional] indicate to LatticeCell the extension of the lattice (via LatticeCell::SetLatticeRange() or LatticeCell::SetNumberOfLatticeElement())

    5. fill the lattice with all the desired Shape_ptr Bi and/or the PinCell Pi (via LatticeCell::FillLattice())

4.8.1 An implicit lattice example (SimpleLattice.cxx & SimpleLattice_serpent.cxx)

This example illustrates implicit lattices, that is to say a lattice filled by a unique universe. A simple tube (a Vessel) is filled by an hexagonal lattice. Each hexagon of the lattice is made of water and filled with a small cylinder of UOx fuel (a pin) inside an Iron cladding. For this example, 2 different implementation are propose to illustrate the only use of shape (like in the 1.d of the beginning of this section) or the use of PinCell (like in 2.a ).

Table 4.1: SimpleLattice example. Implementation using only shapes (left) and using PinCell (right)

//
// define materials (see § 5.1)
//
   Material *H2O= new Material();
   H2O->SetDensity(1.);
   H2O->AddNucleus(1, 1, 2.);
   H2O->AddNucleus(8, 16, 1.);
 
   Material *Iron=new Material();
   Iron->SetDensity(7.87);
   Iron->AddNucleus(26,56,1.);
 
   Material *UOx= new Material();
   UOx->SetDensity(10.2);
   UOx->AddNucleus(92, 235, 0.1);
   UOx->AddNucleus(92, 238, 0.9);
   UOx->AddNucleus(8, 16, 2.);
//
// Geometrical Data Size 
//
   double VesselH=1;
   double VesselR=1;
 
   double HexaH=VesselH;
   double HexaSide=0.2;
 
   double PinR=0.12;
   double CladR=PinR+0.05;
//
//Shapes
//
   //the Vessel : a full Tube
   Shape_ptr Vessel(new Tube(VesselH/2,VesselR));
   //the Exterior
   Shape_ptr Exterior(!Vessel);
   //the lattice generator
   Shape_ptr LatGene(new Hexagon(HexaH/2,HexaSide));
   LatGene->SetUniverse();
 
   // put the lattice generator in the vessel
   LatGene >> Vessel;
 
   // a small cylinder that will contain UOx
   Shape_ptr Pin(new Cylinder(PinR));
   Pin->SetUniverse();
 
   // a bigger one that define fuel clad
   Shape_ptr Clad(new Cylinder(CladR));
   Clad->SetUniverse(Pin->GetUniverse());
 
   //Outside of the Clad
   Shape_ptr OutClad=!Clad;
   Pin>>Clad;
 
//
//Cells
//
   Cell *pin=new Cell(Pin,UOx);
   Cell *clad=new Cell(Clad,Iron);
   Cell *surrounding_clad=new Cell(OutClad,H2O); 
 
   Cell *exterior=new Cell(Exterior,0,0);
   exterior->SetComment("Exterior of the cylinder");
   //vessel
   Cell *vessel=new Cell(Vessel);
   vessel->SetComment("The vessel");
   //the lattice
   //-------------------- Line A --------------------
   LatticeCell *Pavage=new LatticeCell(LatGene);
   //fill the lattice the Pin shape
   Pavage->FillLattice(Pin);
   //-------------------- Line B --------------------

//
// define materials (see § 5.1)
//
   Material *H2O= new Material();
   H2O->SetDensity(1.);
   H2O->AddNucleus(1, 1, 2.);
   H2O->AddNucleus(8, 16, 1.);
  
   Material *Iron=new Material();
   Iron->SetDensity(7.87);
   Iron->AddNucleus(26,56,1.);
 
   Material *UOx= new Material();
   UOx->SetDensity(10.2);
   UOx->AddNucleus(92, 235, 0.1);
   UOx->AddNucleus(92, 238, 0.9);
   UOx->AddNucleus(8, 16, 2.);
//
// Geometrical Data Size 
//
   double VesselH=1;
   double VesselR=1;
 
   double HexaH=VesselH;
   double HexaSide=0.2; 
 
   double PinR=0.12;
   double CladR=PinR+0.05;
//
//Shapes
//
   //the Vessel : a full Tube
   Shape_ptr Vessel(new Tube(VesselH/2,VesselR));
   //the Exterior
   Shape_ptr Exterior(!Vessel);
   //the lattice generator
   Shape_ptr LatGene(new Hexagon(HexaH/2,HexaSide));
   LatGene->SetUniverse();
 
  // put the lattice generator in the vessel
  LatGene >> Vessel; 
 
 
 
 
 
 
 
 
 
 
//
//Cells
//
  PinCell* pinUox=new PinCell;
  pinUox->AddLayer(Uox,PinR);
  pinUox->AddLayer(Iron,CladR);
  pinUox->SetSurroundingMaterial(H2O);
 
  Cell *exterior=new Cell(Exterior,0,0);
  exterior->SetComment("Exterior of the cylinder");
  //vessel
  Cell *vessel=new Cell(Vessel);
  vessel->SetComment("The vessel");
  
  //the lattice
  //-------------------- Line A --------------------
  LatticeCell *Pavage=new LatticeCell(LatGene);
  //fill the lattice with the pincell
  Pavage->FillLattice(pinUox); 
  //-------------------- Line B --------------------

As it can be seen, the version using PinCell is more compact and more easy to understand (Blue part disappear and magenta part is replaced by green one ; the lattice is thus filled by the PinCell). These two implementations give the same result which is shown in Figure 4.2. Note that in the first implementation, outside of the cladding (OutClad Shape) has the same universe than the Clad Shape because it is define as the complementary Shape of Clad (the ! operator). Note also than OutClad and Exterior shapes must be defined before the “put-in” operator (>>) because this latter modify the definition of Vessel or Clad. Few words for Serpent users: the “lat” card of Serpent need an origin, and a pitch. These quantities are deduced from the generator shape ; translating it will thus change “lat” origin. The sizes of the generator shape will determine the pitch (number of values as well as values depending on the dimension (2D or 3D) of the lattice).

pict

Figure 4.2: A simple implicit lattice.

4.8.2 An explicit lattice example with different zones (SimpleLattice2.cxx & SimpleLattice2_serpent.cxx)

Explicit lattices are filled with different universes, allowing to put different structures at some given position. To illustrate this, let us modify the previous example slightly: we want to fill the Vessel with 2 types of hexagons: the first ones are full water hexagons (the reflector) in the outer part of the Vessel and the second ones are the previous hexagons with the small pin rod inside. For shortness, we chose to fill the Vessel with second hexagon type if it is inside a “virtual” cylinder of radius 90cm. Both implementations using Shape only or PinCell has been proposed in Tab. 3.3. Here, a focus is made on the PinCell implementation. The principal changes take place between the lines named “Line A” and “Line B” . We have to replace this part by:

        //-------------------- Line A --------------------
        PinCell* whole=new PinCell;
        whole->SetSurroundingMaterial(H2O);
        LatticeCell *Pavage=new LatticeCell(LatticeGenerator);
        int range=int(VesselR/(1.5 * HexaSide))+1;
        Pavage->SetLatticeRange(-range,range,-range,range); //define a explicit lattice
        for(int i=-range; i<=range; i++)
            for(int j=-range; j<=range; j++)
            {
                int pos[3]={i,j,0}; //the lattice index
                double xt =  HexaSide*sqrt(3.) * (i + j*0.5 ) ;
                double yt = 1.5 * HexaSide * j;
                double X[2]={xt,yt}; //the center of each hexagons
                if(IsHexagonInTube(X,LatticeGenerator,0.9)) //true if the hexa center at X is 
                    Pavage->FillLattice(pinUox,pos);           //inside a tube of radius 90cm
                else
                    Pavage->FillLattice(whole,pos);           //inside a tube of radius 90cm 
           }
        //-------------------- Line B --------------------

The result is shown in Figure 4.3.

pict

Figure 4.3: A lattice with different zones .

Note: if one prefers use only shapes (and not PinCell), one can used the following trick (but it is not implemented for Serpent!!!)

Shape_ptr Whole(new Node(-1));//-1=intersection, here, of nothing
Whole->SetUniverse();

4.8.3 A lattice with more than one simple shape (Stadium.cxx & Stadium_serpent.cxx)

In this example, we consider a hexahedra lattice. Each brick of the lattice will have parts of “stadiums” on each (x,y) sides (see Figure 4.4, this is the base geometry of the MSRE reactor).

//
//materials
//
    Material *Salt=new Material(); //this is the fuel LiF-ThF4/UF4
    Salt->SetDensity(3.1);
    Salt->AddNucleus(3,7,0.7);
    Salt->AddNucleus(9,19,1.9);
    Salt->AddNucleus(90,232,0.15);
    Salt->AddNucleus(92,233,0.15);
    
    Material *Graphite=new Material();
    Graphite->SetDensity(1.86);
    Graphite->AddNucleus(6,0,1);
//
// Data
//
    double VesselH=0.8; // the core H and R
    double VesselR=0.4;
    
    double SquareH=VesselH; //stadium heigth
    double SquareS=0.2;        //stadium straight line length
    
    double BendR=0.02;        //stadium bend curve radius
    double BendH=VesselH;
    double BrickW=0.04;
    double BrickL=0.08;         
//
//Shapes
//
    //the Vessel
    Shape_ptr Vessel(new Tube(VesselH/2,VesselR));
    //the Exterior
    Shape_ptr Exterior(!Vessel);
    //the lattice generator
    Shape_ptr LatticeGenerator(new Brick(SquareS/2,SquareS/2,SquareH/2));
    LatticeGenerator->SetUniverse();
    //put the lattice generator inside the core
    LatticeGenerator>>Vessel;
    
    //Definition of the stadiums (horizontal & vertical)      
    
    Shape_ptr Circle0(new Tube(BendH/2,BendR));
    Circle0->SetUniverse();
    //clone circles for horizontal and vertical stadiums
    Shape_ptr Circle1=Circle0->Clone();
    Shape_ptr Circle2=Circle0->Clone();
    Shape_ptr Circle3=Circle0->Clone();
    //translate circle
    Circle0->Translate(-0.1,+0.04,0);
    Circle1->Translate(-0.1,-0.04,0);
    Circle2->Translate(-0.04,0.1,0);
    Circle3->Translate(+0.04,0.1,0);
    //2 bricks for horizontal and vertical stadiums
    Shape_ptr Brick0(new Brick(BrickW/2,BrickL/2,SquareH/2));
    Brick0->SetUniverse(Circle0->GetUniverse()); //belong to the same universe than Circle0
    Brick0->Translate(-0.1,0,0);
    Shape_ptr Brick1(new Brick(BrickL/2,BrickW/2,SquareH/2));
    Brick1->SetUniverse(Circle0->GetUniverse()); //belong to the same universe than Circle0
    Brick1->Translate(0,0.1,0);
    //2 vertical stadiums
    Shape_ptr Stade0=Circle0 | Brick0 | Circle1; //union of circles and brick
    Shape_ptr Stade2=Stade0->Clone();
    Stade2->Translate(0.2,0,0);
    
    //2 horizontal stadiums
    Shape_ptr Stade1=Circle2 | Brick1 | Circle3;
    Shape_ptr Stade3=Stade1->Clone();
    Stade3->Translate(0,-0.2,0);
    //interior of all stadiums
    Shape_ptr AllStades=Stade0 | Stade1 | Stade2 | Stade3; //to define all stadium exterior
     
//
// Cells
//
    Cell *exterior=new Cell(Exterior,0,0);
    exterior->SetComment("The exterior of the reactor");
    //vessel
    Cell *vessel=new Cell(Vessel);
    vessel->SetComment("The reactor vessel");
    //Lattice
    LatticeCell *Pavage=new LatticeCell(LatticeGenerator);
    Pavage->FillLattice(AllStades);
        
    Cell *stade0=new Cell(Stade0,Salt);
    Cell *stade1=new Cell(Stade1,Salt);
    Cell *stade2=new Cell(Stade2,Salt);
    Cell *stade3=new Cell(Stade3,Salt);
    Cell *exterior_stade=new Cell(!AllStades,Graphite);

pict
Figure 4.4: A lattice with more than one shape.

4.8.4 Lattice of a Lattice (LatticeOfLattice.cxx & LatticeOfLattice_serpent.cxx)

The purpose of this example is to illustrate a hexagonal Lattice of big structural hexagons, themselves filled with a hexagonal lattice of fuel rods (like in a SFR core)

//
// Material
//
    Material *Iron=new Material();
    Iron->SetDensity(7.87);        //density in g/cm3
    Iron->AddNucleus(26,56,1.);
 
    Material *Na=new Material();
    Na->SetTemperature(600);
    Na->SetDensity(0.97);
    Na->AddNucleus(11,23,1);
    Material *Na_cold=Na->Clone(300); //just to see in Serpent all hexagons
    
    Material *Fuel=new Material();
    Fuel->SetTemperature(900);
    Fuel->SetDensity(19.);
    Fuel->AddNucleus(92,238,0.8);
    Fuel->AddNucleus(94,239,0.2);
//
// data size
//
    double VesselH=1.1; //the core size
    double VesselR=1.1;
    
    double FuelHexaSide=0.05; //(small) fuel hexagons side
    double PinRadius=0.02;      //fuel pin radius
    
    double StructHexaSide=8./sqrt(3.)*FuelHexaSide; //(big) structure hexa
    double StructThick=0.01;    //the thickness of the structure (inner
                                // part of this structure is filled with
                                // FuelHexa.
    
//
// Shapes
//
    //the Core
    Shape_ptr Vessel(new Tube(VesselH,VesselR));
    //the Exterior
    Shape_ptr Exterior(!Vessel);
    
    //the outside structure hexagon
    Shape_ptr StructHexa(new Hexagon(VesselH,StructHexaSide+StructThick));
    StructHexa->SetUniverse();
    //the inside structure hexagon    
    Shape_ptr InnerStructHexa(new Hexagon(VesselH,StructHexaSide));
    InnerStructHexa->SetUniverse();
    // put the structure hexagons in the vessel
    StructHexa>>Vessel;
    //
    //a "fuel hexagon"
    //
    Shape_ptr FuelHexa(new Hexagon(VesselH,FuelHexaSide));
    FuelHexa->Rotate(-Pi/2);
    FuelHexa->SetUniverse();
    //put the fuel hexagons in the structure (inner part)
    FuelHexa>>InnerStructHexa;
    
//
// Cells
//
    PinCell *Reflector=new PinCell();             //Reflector=beetwen Core and big
    Reflector->SetSurroundingMaterial(Na_cold);  //hexagons. This is "cold Na" mat
    
    PinCell *FuelHexaBorder=new PinCell();            //Fill the inner part of big
    FuelHexaBorder->SetSurroundingMaterial(Na_cold); //hexagons with cold Na"
    
    PinCell *Carandache=new PinCell();        //This is a fuel pin with
    Carandache->AddLayer(Fuel,PinRadius);     //a cylinder of fuel
    Carandache->SetSurroundingMaterial(Na);     //surrounded by Na
    
    Cell *exterior=new Cell(Exterior,0,0); //out side of the core
    exterior->SetComment("The exterior of the reactor");
    Cell *vessel=new Cell(Vessel);    //the core
    vessel->SetComment("The reactor vessel");
    Cell *innerstructhexa=new Cell(InnerStructHexa);
    innerstructhexa->SetComment("The Inside struct Hexa");
    Cell *exterior_innerstructhexa=new Cell(!InnerStructHexa,Iron); //space beetwen InnerStructHexa
    exterior_innerstructhexa->SetComment("The Big Hexa claddings");//and StructHexa fill with Iron
                                                                    
    LatticeCell *structhexa=new LatticeCell(StructHexa);
    structhexa->SetComment("The Ouside struct of Big Hexa (lattice)");
    int range1=int(VesselR/(1.5 * StructHexaSide)+1);
    structhexa->SetLatticeRange(-range1,range1,-range1,range1);
    double StructHexaWidth=sqrt(3.)*StructHexaSide;
    for(int i=-range1; i<=range1; i++)
        for(int j=-range1; j<=range1; j++)
        {
            int pos[3]={i,j,0};
            double xt =  StructHexaWidth * (i + j*0.5 ) ;
            double yt = 1.5 * StructHexaSide * j;
            
            double X[2]={xt,yt}; 
            if(IsHexagonInTube(X,StructHexa,VesselR))
                structhexa->FillLattice(InnerStructHexa,pos);
            else
                structhexa->FillLattice(Reflector,pos);
        } 
        
    LatticeCell *fuel_hexa=new LatticeCell(FuelHexa);    
    fuel_hexa->SetComment("The fuel hexagonal lattice");
    int range=int(StructHexaSide/(1.5 * FuelHexaSide)+1);
    fuel_hexa->SetLatticeRange(-range,range,-range,range);
    
    double FuelHexaWidth=sqrt(3.)*FuelHexaSide;
    for(int i=-range; i<=range; i++)
        for(int j=-range; j<=range; j++)
        {
            int pos[3]={i,j,0};
            double xt =  1.5 * FuelHexaSide * i ;
            double yt = FuelHexaWidth * ( j + i*0.5 );
            double X[2]={xt,yt};
            if(IsHexagonInHexagon(X,FuelHexa,InnerStructHexa))
                fuel_hexa->FillLattice(Carandache,pos);
            else
                fuel_hexa->FillLattice(FuelHexaBorder,pos);
        }

pict
Figure 4.5: Lattice of a Lattice.

Chapter 5
Materials, Sources, Tallies, ...

5.1 Definition of material

A Material is a set of nuclei with their proportions (default units is %mol), it has a density and a temperature.

For example to define B4C of density 2.52 g∕cm3at 296 K:

Material *B4C=new Material();
B4C->SetDensity(2.52); //default units are g/cm3
B4C->AddNucleus(5,10,0.199); //Boron is composed 19.9% of 510B
B4C->AddNucleus(5,11,0.801); // and 80.1% of 511B
B4C->AddNucleus(6,0,0.25);   // This is natural C (Z=6, A=0) ; for 1 mol of B one has 0.25 mol of C

Just before wrintting Material in MC file, proportions are normalized to 1. The Material::AddNucleus() method has the general form

Material::AddNucleus(Z,A,I,prop);

or

Material::AddNucleus(Z,A,prop);

where Z is an integer corresponding to the proton number of the nucleus, A, an integer for the number of nucleons, I, an integer for the isomeric state and prop, a double for the proportion of the nucleus in the Material. For natural isotopes A=0 ; for ground state I=0 (this is the default value when not specified). See § 5.1.3 for units.

5.1.1 Clone

A Material can be cloned via the Material::Clone(double Temperature) method. There are 2 types of clones depending on:

5.1.2 Mix (Mix.cxx)

Two Materials of different densities can be mixed together to create a third material. This is useful in many situations, e.g., the case of reactor fuels with a given enrichment of fissile isotopes. The principle is the following:

M12=M1->Mix(M2,part,kpMOL);

M12 is a new material and is created from the molar part mixture of materials M1 and M2 (for example Uranium and Plutonium Oxides). M1 and M2 can not be used in the problem, because the Material::Mix method forces them to be fictitious materials by using the Material::ForbidPrint() method: THIS IS DONE IN THE Material::Mix() METHOD. The density of the new material M12, is calculated as

------1----- ρ12 = ω1∕ρ1 + ω2∕ρ2

where ωi and ρi are the mass fraction and density of material i. The temperature of the new material is the same as that of material M1.

5.1.3 Units

Density

Default units (k=constant, d=density) for density is kdGCM3 (g∕cm3); but, to avoid confusion, you must set the density unit in Material::GetDensity() (either kdGCM3 (g∕cm3), kdATCM3 (at/cm3) or kdATBCM (at/barn.cm)). User may give the density in atom/barn.cm by using Material::SetDensity(density, kdATBCM).

NOTE: If the density is not given (i.e. set to the default value 0), it could be automatically found if the proportions are given in atoms/barn.cm or in atoms/cm3 .

Proportions

Default units for proportions are in molar percentage. To give proportion in other units, use

Material::AddNucleus(Z,A,I,prop,UNITS) 

where UNITS stands for (k=constant, p=proportion): kpMOL ("%mol" (default)) or kpMASS ("%mass" ) or kpATBCM ("at/barn.cm") or kpATCM3 ("at/cm3").

5.1.4 Material extension for MC

In MCNP or Serpent, material codes have extensions (like “.60c”). These extensions are very important because they distinguish isotopes at different temperatures, bases, versions,… In practice, if no extension is specified in MCNP file, the first isotope encountered in the xsdir is chosen regardless the temperature imposed in the cell where the material is (but this is not the case in Serpent). In MURE user has different ways to specify an extension.

NOTE:

S(α,β) Treatment and MT card

When the special treatment S(α,β) is desired for low energy, MCNP MT card arguments or Serpent moder and therm cards have to be specified. Consider the water example:

H2O->AddNucleus(1,1,2.,kpMOL,"H2O"); 
H2O->AddNucleus(8,16,1.); 

5.1.5 Pseudo Materials

PseudoMaterials are particularly useful for thermal-hydraulics problems. Normally, the material temperature dictates which extension (e.g. “.60c”) for cross-sections will be used from the evaluated data libraries. The default behavior in MURE is that the nearest available temperature is used, the maximum tolerable difference between this temperature and the desired material temperature being defined in the TemperatureMap class. However, there are occasions where it is useful to perform a stochastic interpolation at the desired temperature by adding two identical nuclei at different temperatures in the material in varying proportions. These types of materials are called PseudoMaterials, because, in principle, they behave as a material at a precise temperature, T, even though the cross-sections in the evaluated data libraries at this precise temperature do not exist. To declare a Material (Mat) to be a PseudoMaterial requires the following simple line:

Mat->SetPseudoMaterial(); 

From this point onward, MURE will automatically create two identical nuclei with different extensions/temperatures and in the correct proportions for, every nucleus which is added to this material. No other commands are necessary. However, in the case of an evolving material (a reactor fuel rod for example), this can cause the addition of many, many extra nuclei to the material, thus slowing down MCNP calculations by perhaps as much as 20%. In the case that the user only wants to declare certain nuclei to be PseudoNuclei (for example, the most important ones, such as the fissile and fertile nuclei) the following kind of lines can be added after the creation of a Material:

UO2Fuel->AddNucleus(92,235,0.04); 
UO2Fuel->AddNucleus(92,238,0.96); 
UO2Fuel->AddNucleus(8,16,2.0);
UO2Fuel->GetNucleus(92,235)->SetPseudoNucleus();
UO2Fuel->GetNucleus(92,238)->SetPseudoNucleus();

Now, only the selected nuclei will be treated as pseudo nuclei and automatically inserted twice into the MCNP file at the two nearest data base temperatures in the correct proportions. The proportions, ω1 and ω2 , for the interpolation are calculated in the following way:

√ -- √-- --T-−--T2- ω2 = √T2-− √T1-

ω1 = 1 − ω2

where T1 and T2 are the nearest two temperatures in the evaluated data base files straddling the desired temperature T.

5.1.6 MC material Printing

This section concerns only evolving Materials. When evolution is required, at the beginning, many nuclei are present in the material but in a very small proportion. Each time a neutron enters in a cell, MC code tries to find the nucleus on which the potential reaction could happen. Thus, the more complex is a material, the longer is the MC run. The idea is to suppress only for the MC print, the nuclei which do not contribute “significantly” to neutron transport to save CPU time. The criterion is based on the evaluation of total cross-section σT and nucleus proportion p. If a nucleus of a material has a p × σT < 𝜀, it is not printed. The value of 𝜀 is by default set to 0.1% × 1mb = 10−6. This can be changed via MURE::SetMCNPNucleusThreshold(). To be noticed that perturbative tallies of all evolving nuclei are written (in order to evaluate the cross-section...). The effect of this procedure is thus to accelerate MC runs at the beginning of the evolution where time steps between MC are generally smaller, and add more and more nuclei as far as the evolution is going on (with larger time steps between MC runs).

5.1.7 Examples

Please take a look at Material_test.cxx and Material_test2.cxx and the corresponding Serpent version for the second one. The use of Nucleus::SetXSExtension or Material::SetDefaultXSExtension supposed that the user provide a correct xsdir file with all necessary information.

5.1.8 Automatic extension finding

Here is a short description of the automatic nucleus extension finding.

Making the BaseSummary.dat file

In order to work, a specific file, BaseSummary.dat, must be present in the data directory. The aim of this file is to facilitate the search for a given Z,A,I and related xsdir/xsdata information (see also § 1.4.2 and 1.4.3).

Extracting the MC nucleus code from the BaseSummary.dat file

Suppose the desired ZAI exists in the BaseSummary.dat file. Then the choice of the desired nucleus is based on 2 criteria: the first one is related to the preference in a nuclear base (ENDB, JEFF, JENDL) and the version of this base, and the second criterion is related to the nucleus temperature. A score is calculated for each criterion, and the best score is chosen. Two classes are in charge of these criteria: BasePriority and TemperatureMap.

5.1.9 Automatic XSDIR construction

Using the extension finding, MURE can build the XSDIR/XSDATA file that is used by MC code to get the cross-sections used in the problem. To use this option, you just have to set :

gMURE->SetAutoXSDIR();

When getting the code for MC file, MURE will then save the best line it finds in BaseSummary.dat. Therefore, right before building the MC file (the XSDIR is built at the first step of an evolution), MURE will build the XSDIR/XSDATA file by using all these best lines it saved.

IMPORTANT : If an existing XSDIR/XSDATA is already present in the MC run directory, then MURE will NOT build a new one ; the existing one will be used by MC code. Thus either the xsdir/xsdata is correct, either you have to remove it.

5.2 Particle Source

Because the difference between particle source definition in MCNP and Serpent is great, MURE source definition uses 2 different classes, MCNPSource and SerpentSource to avoid confusion for the user ; they are sibling of the MCSource class.

The MCNPSource SerpentSource classes enable the definition of particle sources (up to now only neutron sources may be defined for Serpent). Nevertheless, some source definition methods are common to both classes. Two type of source can be defined, source for criticallity calculation (called kcode mode) and fix sources for non fissile or subcritical system only. For Serpent, source possibilities are much more simpler than for MCNP ; in MURE, the source implementation for Serpent is also much less developed. Mainly position and energy of punctual monoenergetic isotrope source can be defined. Criticallity is defined also by the “SerpentSource::SetKcode” method. In the following section, the MCNPSource is used because of its greater complexity, but looking to SerpentSource methods and their analog counterpart in MCNPSource is self understandable.

The source may be defined with MCNP card or as an external source. In the latter case (MCNPSource::SetExtern()), external source will be read in a special file. This is valid in subcritical mode as well as in critical (KCODE) mode (not true for Serpent). One chooses to run in subcritical mode (default) or in KCODE mode via the MCNPSource::SetKcode() method.

5.2.1 Setting Source for MURE

Once an MCSource has been defined, one has to give this source to MURE class:

MCNPSource *mysource=new MCNPSource(1000); //1000 source particles 
gMURE->SetSource(mysource);

5.2.2 Criticallity calculation: Kcode mode

To be noticed that, for Kcode, the number of source particles defined in the constructor of MCSource is in fact the number of neutrons per cycle.

To initiate a KCODE, a SDEF (or KSCR) card is needed (if no external source is given!). By default in MCNP, the position of this initial source is (0,0,0). If you want to change the source position, use the MCNPSource::SetPosition(). In Serpent the default positions are chosen in cells containing fissile materials.

MCNPSource *mysource=new MCNPSource(1000); //1000 source particles/cycle
//Kcode with first expected keff=1, 10 inactive cycles and 80 active ones
mysource->SetKcode(80,10,1); 
gMURE->SetSource(mysource);

If one want to use a source from a previous KCODE:

mysource->UsePreviousRunSource(“olds”);

where the source file olds is the result of a previous run. In case of evolution, if you do not specify any argument, the source will be calculated as specified by inactive cycles, and then for next MCNP runs, the source is taken from the previous cycle.

5.2.3 More elaborated sources (MCNP only)

The MCNPSource gives more details on MCNPSource class possibilities. Two simple examples are presented:

First

MCNPSource *source=new MCNPSource(1000); //1000 particles per cycle
source->SetKcode(500,20,1); //critical mode (500 active cycles, 2O inactive, supposed keff=1)
source->SetKSRC(); //use fission spectrum for neutron energy
source->AddPosition(0,0,0); //add one position for the source at (0,0,0)
source->AddPosition(0.1,0.1,0.1); //add an other one at (10cm,10cm,10cm)
source->AddPosition(0,0,0.2); //add an other one at (0,0,20cm)
source->AddPosition(0,0,-0.2); //add an other one at (0,0,-20cm)

and the second

MCNPSource *source=new MCNPSource(); //SDEF source with 10000 particles
source->SetEnergie("D1");  //energy is define by the source probability 1 
source->SetDistribution("RAD=D2 EXT=D3"); // radial and axial extension are 
                                          // defined by source information 1 and 2
source->AddBias("SP1 -5 2.0");  //p(E)=E*exp(E/2.MeV)
source->AddBias("SI2 1. 2.");   // radial extension from r=1cm to r=2cm
source->AddBias("SI3 -10. 10.");// axial extension from z=-10cm to z=10cm

5.2.3.1 Source define via Spectrum class (MCNP only)

One can define energy distribution source using the Spectrum class and the MCNPSource::UseThisEnergyDistribution() method. Any call to that method will dump in MCNPSource_YourEnergyDistribution.dat, the ASCII spectrum of this source (for user own use) as well as the total number of particles emitted (for normalization purpose).

To be noticed: this method will change the source particle type according to the Spectrum type (neutrons for NeutronSpectrum, γ for GammaSpectrum, …).

This method only change energy distribution of the source (using MCNP distribution source cards (D, SI, …) numbered by default from 800).

5.2.3.2 Tube Source (MCNP only)

It is possible to define volumic cylindrical sources with the help of the Tube class. Two methods are available : MCNPSource::AddTubeSource(Tube *tube, string energy) and MCNPSource::AddTubeSource(Tube *tube, Spectrum *spectrum). In both cases, the tube is used to define the shape of the source. The first one allows to use either mono-energetic sources or using source distributions “D” card. In the latter case, the user need to provide himself the distribution description via the MCNPSource::AddBias() method (see TubeSource.cxx example). The method using Spectrum class argument use Spectrum to define both particle type of the source (neutrons for NeutronSpectrum, γ for GammaSpectrum, …) and energy spectrum (see TubeSource2.cxx example).

Theses 2 methods use internally source distribution cards (D, SP, SI, …). The starting number of these distribution is 900. Thus don’t use number above 900 to define your own source distribution and biasing.

Of course you can add more than one Tube source ; but it is only possible to add more Tube sources with same calling method (either MCNPSource::AddTubeSource(Tube *tube, string energy) or MCNPSource::AddTubeSource(Tube *tube, Spectrum *spectrum). One cannot mix them.

5.2.3.3 Define MCNP Source from the result of an evolution (MCNP only)

If one want to use the energy distribution of a spent fuel which has cooled for a time t, use MureGui and dump the source in a MURE format (“Save Data” button, then select the “MURE INPUT Spectrum” radio widget). The generated file can be inserted in a MURE “cxx” code (c.f. MureGui section 8.3.6).

5.3 Tally class

The Tally class enables tallies for MCNP/Serpent to be defined.

WARNING : NEVER USE : “Tally *f=new Tally[N]” BUT USE INSTEAD “vector<Tally*> f(N) ; for (int i=0 ; i<N;i++) f[i]=new Tally;”. The reason is that when tallies are deleted, one need to use either “delete []” or “delete” ; in MURE, the second case is chosen which is incompatible with the former declaration.

MCNP tallies or equivalent Serpent detectors are used to store quantities from a run ; Serpent detectors and MCNP tallies have lots of common points but the MCNP version has more possibilities (Group tallies, tally normalization, ...). The implementation of tallies in MURE is (at least for historical reasons) much more developed for MCNP. Nevertheless, surface tallies (type 1 and 2 in MCNP) are not fully implemented (Surface calculation, ...) even if one can use them with special precaution/verification.

There are 2 main tally types: Surface tallies (type 1 and 2 in MCNP, same convention is kept for Serpent) and Cell tallies (4,...). At present, only Cell tally of type 1 (current through a surface), 2 (flux through a surface), 4 (flux in cell), 6 (energy deposition in a cell) and 7 (fission energy deposition in a cell) are implemented in MURE. Surface tallies have not been tested at all and stochastic surface calculation is not implemented in MURE.

A Tally consists of

Of course Cell tally bin types could not be mixed with surface tally bin type. Tally bins are added with the Tally::Add() method (for cell/surface bins) whereas energy and/or time bins are added with Tally::AddEnergy() or Tally::AddTime().

Tally multiplicator may be added with Tally::AddMultiplicator() to obtain cross sections, reaction rates and so on.

Example of Tallies:

Tally *t1=new Tally(); //default tally type is 4: flux in cell
t1->Add(cell_1);  
t1->Add(cell_2);
LatticeBin *lpb=new LatticeBin(rod); //rod is a cell associated to an element of a Lattice
lpb1->AddContainer(core); //core is a cell fills by "rods"
t1->Add(lpb);
t1->AddEnergy(1e-2,1.e7); //log energy binning from 0.01 eV to 10 MeV with 10 bin/decade (default)
t1->AddMultiplicator(mat_1,new Reaction(102)); //capture reaction rate
t1->AddMultiplicator(mat_1,-6);//fission reaction rate

 

When in LatticeBins, you can precise the positions of the cell in the lattice, use a string which gives the position in the MCNP way, that is :

string pos="[0 0 0]";
LatticeBin *lpb=new LatticeBin(rod); //rod is a cell associated to an element of a Lattice
lpb->AddContainer(core,pos); //core is a cell fills by "rod" (e.g. a lattice)

This means that you are interested by the flux in the cell rod which is in the cell core at the lattice position [0,0,0]. The result in MCNP file will be (given that the number of cell rod is 1 and cell core is 2):

F4:N (1<2 [0 0 0])

Finally, you can use the universe shorthand (u=1 for instance), simply by giving the universe number, either directly if you know it or by using Shape::GetUniverse() :

Tally *t1=new Tally();  
t1->Add(MainShape->GetUniverse());

Please note that, in MCNP, the universe shorthand doesn’t work for lattice elements filled with the universe number, that is to say, when a lattice is filled by position, for instance when in MURE we use Cell::FillLattice(Universe_Number, Position), the universe is not detected by MCNP. (see example Tally_Test.cxx)

As this shorthand is replaced by all the cells filled with the universe, we strongly recommend that you use it only if you are perfectly sure of the number of cells involved in the tally defined in that way.

Each tally must have a definite volume (or surface). The volume is set via the Tally::SetBinVolume() or for a given bin, via the TallyBin::SetBinVolume() or TallyBin::SetVolume(). This may be done manually or automatically: if no bin volume has been provided, the volume is automatically calculated by a stochastic method (see also the note C.1 concerning the volume calculation). For lattices and universe shorthand, as multiple volumes are needed for each bin, the SetBinVolume() method precises a partial number.

Tally *t1=new Tally(4);
t1->Add(Cell1);
t1->SetBinVolume(0,1.e-6,0); // Set the volume of Cell 1 to 1 cm3
LatticeBin *lpb2=new LatticeBin(rod1);
lpb2->AddContainer(rod2);
lpb2->AddContainer(core);
Tally *t2=new Tally(4);
t2->Add(lpb2);               // A Lattice Bin giving somthing like (1 2 < 3)
t2->SetBinVolume(0,1.e-6,0); // Set the volume of Cell 1 in cell 3
t2->SetBinVolume(0,1.e-6,1); // Set the volume of Cell 2 in cell 3

The partial number follows the order of the output. (to know exactly the order for lattice bin volumes, see MCNP doc or example Tally_test.cxx).

5.3.1 Fluence to dose conversion (MCNP Only)

The Tally::AddFluenceToDoseConversion() method converts fluence (only for F4 type tally) to dose rate (rem/h) using fluence to dose conversion factors. It can be use only for photons or neutrons. The conversion factor are store in /path_to_mure/MURE/data/GammaFluenceToDose.dat and NeutronFluenceToDose.dat (they are taken from [23]). Interpolation is done for energies and for factors. This interpolation is done by default in log way. The 2 boolean parameters of Tally::AddFluenceToDoseConversion() allows to specify linear interpolation for energy (the 1st one) and conversion factors (the 2nd one).

Tally *MyTally=new Tally(4,ParticleType);//fluence to dose works only for ParticleType=”P” or “N” 
MyTally->Add(....);
MyTally->AddEnergy(.....);
MyTally->AddFluenceToDoseConversion();

Chapter 6
Nuclei Tree

6.1 Nuclei Trees: general considerations

To evolve a given material, it is essential to know in advance the ensemble of nuclei which it can produce via successive nuclear reactions and decays. This information comes from the Nuclei Tree object, which is defined by the NucleiTree class, and contains the linking information between all the nuclei that exist in the chart of nuclides. Nuclei can be transformed into other nuclei by nuclear decays (β, α, electron capture, etc.) and nuclear reactions ((n,γ), (n,2n), fission, etc.). In general a particular nucleus has daughters (those nuclei it can transform into) and parents (those nuclei that can produce it). To calculate how much of a particular nucleus is produced during the evolution it is essential to know the ways in which it can be produced and the ways it can be destroyed.

We cannot only define the total global nuclear tree (the connections between all 3834 nuclei), but also the local or particular tree for a given nucleus. The particular tree for a certain nucleus is, in general, a subset of nuclei and their linking from the global tree. To extract an individual tree from the global tree, the NucleiTree::ExtractZAI method is used. The local tree is explored recursively by starting at the given nucleus (see figure 6.1) and then finding each successive reaction and decay daughters.

pict

Figure 6.1: The chart of nuclides used in MURE containing 3834 Nuclei. Nuclei colored in red are those present in the Nuclear Tree of 232Th. To build this tree, simple rules are used for determining which reactions are possible: if Z < 50, (n,γ) reactions are allowed, if Z < 90, (n,γ) and (n,2n) reactions are allowed, and if Z ≥ 90, (n,γ) and (n,2n) and fission reactions are allowed.

6.1.1 How the tree information is stored

pict
Figure 6.2: The general case of a ZAI object X in a tree with multiple reaction/decay parents and daughters.

6.1.2 Cutting or simplifying the tree

Isotopes in the chart of nuclides have half-lives ranging over ~40 orders of magnitude. In a reactor core neutron fluxes are much lower (∼ 1015 neutron/cm2/s) than in supernova explosions (up to ∼ 1024 neutrons/cm2/s) and thus decay rates are usually much higher than reactions rates for most short-lived nuclei. Furthermore, the reactor fuel composition only changes significantly after some years, so short-lived isotopes play no important role in the evolution. We can therefore assume that nuclei with short half-lives decay instantaneously which allows us to simplify the tree and reduce the number of nuclei in the calculation and hence the calculation time. To simplify a nuclear tree requires these short-lived nuclei to be cut out of the tree and their reaction/decay parents and daughters re-linked into the tree correctly. There are two types of simplification for a particular tree branch, depending on whether the nucleus X in question is produced via decay (see figure 6.3) or nuclear reaction (see figure 6.4) .

pict

Figure 6.3: Tree simplification if nucleus X is produced via radioactive decay

Note that this means that in some cases a single type of reaction can produce multiple daughters with different probabilities (branching ratios), since the reaction daughter had more than one decay daughter itself, and was cut out of the tree because its half-life was too short (see figure 6.4). After the simplification process, these reaction daughters will have weights which are proportional to the decay branching ratios of the nucleus that was removed from the tree.

pict

Figure 6.4: Tree simplification if nucleus X is produced via a reaction. The decay daughters of X become reaction daughters of A.

For a particular tree many nuclei towards the neutron drip-line end up getting cut out due to their progressively shorter half-lives. In figure 6.5 we can see the nuclei explored and then removed from the tree right up to the edge of the nuclear chart.

pict

Figure 6.5: The simplified nuclear tree for 232Th. Nuclei in yellow are those which were cut out during the tree exploration process because their half-lives were too short (less than 1 hour). Nuclei in red remain in the tree. Reactions follow the same rules that those defined in fig. 6.1.

6.1.3 The difference/similarity of material trees

Two nuclei X, and Y will have identical trees if there is a pathway of decays and reactions that produces Y from X and also a pathway of decays and/or reactions that produces X from Y. Reactions (n,γ) and (n,2n) move up/down in mass, and decays only move down in mass. Once the mass is sufficiently low that (n,2n) reactions are not permitted, then decays will only take nuclei as far as the valley of stability and a stable nucleus are reached. At this point the tree branch ends. However, successive neutron captures are always possible, therefore fuels containing lighter elements, such as the Oxygen in Uranium Oxide (UOX) fuel, will have a different tree from a pure fuel’s one such as the 232Th (see figure 6.6).

pict

Figure 6.6: The UOX tree for T1∕2 ≥ 1 hour. Reactions follow the same rules that those defined in fig. 6.1.

6.1.4 How many nuclei need to be considered in the evolution?

In all the previous chart figures, the reactions taken into account are: if Z < 50, (n,γ), if Z < 90, (n,γ) and (n,2n), and if Z ≥ 90, (n,γ) and (n,2n) and fission reactions. In reality, of course, nuclear cross sections are only available for a subset of all nuclei. Figures 6.7 and 6.8 show the nuclei with available reaction data found in the default MCNP4 cross-section (from ENDFB base). All the available reactions are plotted in figure 6.9.

pict

Figure 6.7: 232Th tree with MCNP default cross-sections and T1∕2 ≥ 1 hour.

pict

Figure 6.8: The UOX tree with MCNP default cross-sections for T1∕2 ≥ 1 hour.

pict

Figure 6.9: MCNP default available reaction data are shown in light blue.

For a realistic half-life cut, the actinide nuclear tree still contains several hundred nuclei. However, simplification has reduced considerably the complexity of the problem(see figure  6.10)1.

pict

Figure 6.10: Number of nuclei in the 232Th tree as a function of the minimum half life allowed, T1∕2min. The blue curve shows the number of nuclei in the tree if the simple rules for allowed reactions are used (see figure. 6.1). The black curve shows the number of nuclei in the tree if allowed reactions are taken from the MCNP default reaction data base. The green curve corresponds to the MURE data base built at IPNO/LPSC (not part of the MURE package). The red dot shows the number of nuclei in the tree for the default cut (1 hour).

6.1.5 Fission Products

Fission products are built according to the fission yields available in nuclear data libraries (default library is ENDF B 6.8). Fission yields are read in the FPyield.bin with the help of FPavailable.dat files (in the data directory). These files have been generated using SMURE/source/utils/fp/GenerateFPYield2 and can be regenerated using any ENDF fission product yields of ENDFB.X, JEFF.X, JENDL.X, etc libraries.

To be noticed, that for a given heavy nuclei, if the FP yields are not available, the “closest” Z,A,I heavy nuclei FP yields are used (which is of course wrong but probably better than nothing). For example, spontaneous fission of 244Pu uses the fission product distribution of the induced fission of 243Am. If some FP of the libraries are not included in the decay data file (chart.jeff3.1.1), then this fission product is ignored or, for isomeric states, the yield is given to first state below. Up to January 2024, the only case found is the production of the second isomeric state of 116Ag found in JEND4 and JENDL5 which the yield is added to the one of the first isomeric state of this nucleus.

They are the major contributor to nuclei increased in the reaction tree. The more FP in nuclear libraries, the longer the MCNP/Serpent execution time will be due to the increase in tally number. Thus, in order to obtain a first approximation, it may be convenient to run only a selection of fission products that may be considered as “fundamental”. We provide a method, NucleiTree::FissionProductSelection() to choose the most important FP. All other FP will not evolve (no resolution of Bateman equations) but instead are added in the real materials of the MC input file in order to have a good flux description. A default selection has been entered with the nuclei of Table 6.1 ; nevertheless the user can provide a file with its own selection. To use this selection you must use gMURE->KeepOnlyFissionProductSelection() ; to give your own selection, use gMURE->KeepOnlyFissionProductSelection(“MyFPSelection.dat”) where the file MyFPSelection.dat is an ASCII file with each FP you prefer to keep for evolution : for example

Z1 A1
Z2 A2
...

Table 6.1: Default 58 fission product selection in NucleiTree::FissionProductSelection().
36-Kr- 83 46-Pd-107 54-Xe-134 60-Nd-145 62-Sm-152 66-Dy-160
40-Zr- 93 46-Pd-108 54-Xe-135 61-Pm-147 63-Eu-151 66-Dy-161
42-Mo- 95 47-Ag-109 54-Xe-136 61-Pm-148 63-Eu-152 66-Dy-162
43-Tc- 99 48-Cd-113 55-Cs-133 61-Pm-149 63-Eu-153 66-Dy-163
44-Ru-101 49-In-115 55-Cs-134 62-Sm-147 63-Eu-154 66-Dy-164
44-Ru-103 51-Sb-125 55-Cs-137 62-Sm-148 63-Eu-155 67-Ho-165
44-Ru-106 52-Te-127 55-Cs-135 62-Sm-149 64-Gd-155 68-Er-166
45-Rh-103 53-I -127 56-Ba-138 62-Sm-150 64-Gd-156 68-Er-167
45-Rh-105 53-I -135 59-Pr-141 62-Sm-151 64-Gd-157
46-Pd-105 54-Xe-131 60-Nd-143 64-Gd-154 64-Gd-158

6.1.6 Isomer Production from neutron reactions

TO BE NOTICED: Before July 2012, MURE only takes into account the case of 241Am. In the previous treatment, all ground state production was replaced by 242Cm ; this was of course wrong and lead to an over estimation of 242Cm production and an under estimation of 242Pu. This has been corrected since July 2012.

In some cases, neutron reactions like (n,γ), (n,2n) or other reactions can lead to either ground state or isomeric states of a nucleus. From all the many isomeric levels only the long-living are of interest, because for neutron physics application the isomeric state behaves like a new isotope with different half-live, reaction and decay channels. For example, let’s take the case of 241Am. Neutron capture on this nucleus leads to 242Am(T1∕2 ∼ 16h) and 242mAm(141y). Ground state decays in 242Pu(17.3%) by electronic capture and in 242Cm(82.7%) by β− decay. A file (IsomerProduction.dat) is provided to give the branching ratio (for thermal and fast spectrum) to produce its ground and isomeric states of a large set of nuclei. The format of this file is

ZtAtIt ZpApIp ReactionCode BRthBRf halflife Cutflag

for any isomeric state branching of a nucleus, where Z, A, I are proton number, mass number and isomeric state number, respectively. Subscripts t and p indicate target and product nuclei. The ReactionCode denotes the ENDF MT reaction number and BRth and BRf are the branching ratios averaged over a thermal (th) and fast (f ) neutron spectrum. 6.2 shows as an example the production (index p) of Am-242m from target (index t) isotopes 241Am, 242Am, 243Am, 244Am via the capture (102), inelastic scattering (4), n2n (16), n3n (17) neutron reactions, respectively. The branching ratio is given for a thermal LWR neutron spectrum (index th) and a fast spectrum (index f ).

Table 6.2: Isomer production of 242mAm.
Zt At It Zp Ap Ip ReactionCode BRth BRf
95 241 0 95 242 1 102 0.13 0.15
95 242 0 95 242 1 4 0.09 0.06
95 243 0 95 242 1 16 0.26 0.26
95 244 0 95 242 1 17 0.24 0.25

The “halflife” parameter is the half life of the produced ground state (T1∕2 ∼ 16h for 242Am) followed by its unit (here “h”) ; this parameter in not used, it is just an indication that can help the user to decide if the last parameter (CutFlag) is a “X” or a “M”. “M” means that MURE decide what to do according to MURE::GetShortestHalfLife() method. “X” forced the replacement of the ground state (e.g. 242Am) by its decay daughters. To be noticed, if “M” is chosen, wrong results can be obtained in 242Pu and 242Cm for example, if no cross section is available for the ground state (242Am). Currently, the mechanism of “X” and “M” is disabled and MURE handles everything according to MURE::GetShortestHalfLife().

Branching Ratio To use the branching data the corresponding reaction type must be available in the Ace Monte-Carlo library. This might be not the case for some total reactions like inelastic scattering which is the sum of all inelastic levels and continuum stored as MT 51–91. To generate this table, the nuclear data processing code PREPRO [24] is used to generate one group activation cross sections for any isotope of the TENDL 2014 [25] nuclear data evaluation. In PREPRO, the weighting spectrum options of the NJOY [26] code for spectrum type 5 (EPRI CELL) LWR and type 8 fast neutron spectrum were applied. In total, 30562 reactions of almost 3000 isotopes leading to isomeric states are generated. In the course of the depletion calculation MURE requests reaction rates from the MCNP run and applies the branching ratios depending on the target and production isotope to fill the transmutation matrix.

As another source for branching ratios for isomer production the nuclear data service at NNDC-ENDF maybe used, in the “advanced” or “extended” retrieval tab, by choosing “n,g” in the Reaction line and “MRNP” (Multiplicities for production of radioactive elements, MF=9) in the Quantity line. The branching ratio given in MURE/data/IsomerProduction.dat are effective branching ratio computed on a typical thermal and fast spectrum as explained in [27]:

∫∞ -0-B∫R-(E)σreact(E)ϕ(E)dE- BReff = ∞0 σreact(E)ϕ(E)dE

The NNDC-ENDF site gives either the branching ratio to the ground state or to the isomeric state(s) ; in MURE only one branching ratio maybe given (e.g. to the first state), the ground state ratio is then added. However, sometimes more than one branching ratio is needed. In that case, all ratios and reactions have to be added to the file. 2.

6.2 Nuclei Tree: the implementation

6.2.1 Important files

All provided data files are located in MURE/data. One can use either, gMURE->SetDATADIR(path) or the environment variable DATADIR (setenv DATADIR path in csh).

6.2.2 Reaction Auto-detection

Reaction auto-detection is based on information found in the ACE binary4 5 format files. Nuclei, for which reaction cross-section database information is available, can be found. The auto-detection process relies on the ReactionList class.

How does it work?

PLEASE NOTE: if the BaseSummary.dat file is modified (add new nucleus, or new base) or if the threshold used to find significant reactions is changed, THE ReactionList DIRECTORY MUST BE REMOVED. Otherwise, your modifications will not be taken into account.

What are significant reactions?

A quantity Ξi is calculated as

∫ σ dlog E Ξi =-∫-i----- dlogE

where i stands for reaction i. Then if and only if Ξi ≥ σmin where σmin is a minimum threshold, the reaction is take into account ; The default threshold has been set to 0.01 barn6. This can be modified by MURE::SetReactionThreshold() called before any MURE::SetReactionListInitMethod() or Material definition. In order to save time, the procedure to allow or disable a reaction for a given nucleus is performed only the first time an evolution for that nucleus is asked ; then the result for that nucleus (i.e., a ReactionList object) is stored in a file (one per Z) in a specific local directory named ReactionList.

NOTE Any change of the MURE::SetReactionListInitMethod() or MURE::SetReactionThreshold() after a first run will be taken into account only if the nucleus has not been already processed: thus you must destroy the ReactionList directory in order to have these changes taken into account.

6.2.3 The recursion depth cut

It is possible to put a limit on the recursion depth of the reactions when calling the NucleiTree::ExtractZAI method. For example prohibiting 30 or more successive neutron captures limits the nuclei included in the tree to the top of the actinide region (e.g Fermium, Californium, etc.) where experimental cross-section information does not exist, and many successive neutron captures are the only way in which these nuclei can be produced. The recursion depth parameter should be used with care. Figure 6.13 shows the effect of the recursion depth. The default value is 10000...

pict

Figure 6.13: Effect from limiting the recursion depth for successive reactions in the 232Th tree, where nuclei with T1∕2 < 1 hour have been removed.

Chapter 7
Evolution in MURE

7.1 Preliminary Remark

Before doing any evolution you must provide a system (geometry based on Cells, all Materials, running MC code mode (critical or subcritical), ...). This system definition is normally provided by the MURE geometry input file (a C++ code). But if you have a very large MCNP file already defined, it is possible to perform an evolution using this file without redefining the whole MCNP geometry in MURE style. The way of doing this is explained in section 7.5. This has not been yet implemented for Serpent.

7.2 General Considerations

Simulating time evolution involves solving the well-known Bateman equations for every given nucleus in every given cell for which evolution is desired, and carrying out this integration over a discrete number of time steps. Simulating the fuel burn-up may also involve adjusting the neutron flux value at every step of the integration process so as to model constant reactor power, adjusting poison concentration to keep a given value for keff, ...

Implementation of evolution related methods are handle via the EvolutionSolver class. An object of that class is automatically built in MURE constructor, thus on can access to it via gMURE->GetEvolutionSolver().

The general scheme for the evolution is shown in figure 7.1:

pict

Figure 7.1: Evolution principle.

At the present time:

The results of an evolution consist of a number of files stored in an “evolution” directory ; the name of this directory should be defined by the user via the MURE::SetMCRunDirectory function. Two general sorts of files can be found in this directory : On the one hand, there are the MC input and output files (“o”, “m”, “s” and so on for MCNP and “_res.m”, “_det0.m”, ... for Serpent) and on the other, the MURE-generated summary files (DATA_i, BDATA_i, ...) that can be read by a ROOT (see ROOT home page) C++ interface: MureGui.

7.2.1 Time discretization

There are 2 or 3 levels of time discretization:

7.2.2 The Bateman’s Equation

Here is the general form of the first-order differential equation that has to be solved for every nucleus in every evolving cell.

∑ ∑ ∑ ∂Ni-= − λiNi + λj→jiNj + NjRj →iσr,jϕ − Ni σr,iϕ ∂t j j,r r ◟-------◝◜-------◞ ◟-----------◝◜-----------◞ Decay Reaction

In the decay part, the meaning of j → i is that the sum of j is done only for nuclei j that decay to i ; in the production of the reaction part, the branching ratio is 0 if the reaction r on j does not produced i, 1 in most of production reaction, except for fission where Rj→i is the fission yields of nucleus i associated to the fission of j, or when reaction produced an isomer and a fundamental state of a nucleus, then it is the production branching ratio of the wanted state (in a evolution, isomer and fundamental are considered as 2 different nuclei).

The writing and solving of these equations are systematically carried out by using internal methods of the EvolutiveSystem class. Only the general approach will be specified here. This approach involves the definition of what is called the evolution matrix E for every evolving cell. The user does not need to know in detail how it works and it is only being described here for the record.

7.2.2.1 The evolution matrix

Internally, every evolving cell knows which nuclei it contains (initially) or will eventually contain due to the evolution of its initial composition, thanks to the prior construction of the nuclei tree (c.f. NucleiTree in chapter 6). The approach taken in MURE is to assign unique indexes to every nucleus in the cell’s composition so as to build its evolution matrix E. It is a square matrix of n dimensions, n being the total number of nuclei in the composition of the cell.

After completion of the matrix, the ith nucleus of a given cell will evolve according to its own Bateman’s equation which will then be simply written as the sum over all nuclei:

∂Ni-= ∑ E × N ∂t j ij j

where Eij is the corresponding element of matrix E and represents the contribution of the jth nucleus in the cell’s composition to the rate of change in time of the ith nucleus.

The actual writing of the matrix follows the natural structure of Bateman’s equations and is done in two separate parts : first by taking into account all parent/daughter decay relations between nuclei considering proper decay branching ratios. After this, the parent/reaction daughters term is added, where, in the case of fission, fission yields are properly considered.

The evolution matrix is rebuilt in this way at every CRAM/Runge-Kutta sub-steps (c.f : 7.2.1 Time Discretization) for all evolving cells. The solving for the material’s composition new proportion vector −→N = {N 0, N1 ,…,Nn} is done automatically using either a CRAM or a fourth-order Runge-Kutta-type integration method.

7.2.2.2 Multi-Threading parallelization

Since gcc version 4.1, OpenMp multi-threading is available with gcc compilation. Depending on the distribution, libgomp library (omp.h and libgomp) need to be installed.

The multi-threading behavior is controlled by the OMP_NUM_THREADS environment variable. Setting it to 1 is equivalent to the single processor mode. Setting it to the maximum number of cores available, will run on all that cores (which is the default when this variable is not set). This variable only controls the behavior of the evolution NOT the MC run itself (which is still controlled by the MURE::SetOMP() method).

7.2.2.3 Cooling period

Whereas in general the evolution is useful when the reactor runs at its given power, it is also possible to introduce some cooling period, where the power (and thus the flux) is set to 0. This is done by using an array of booleans of same size than the time vector at which MC runs are normally performed. Each time this array contains a “true”, the flux is imposed to 0 for the considered step, else the evolution is a standard one. The array is passed to MURE via gMURE->GetEvolutionSolver()->SetCooling() method.

7.2.3 Solution of the Bateman equation

The numerical solution of the first order differential Bateman equations is subject to research for a long time. Since the half-lives vary from some seconds to million of years the stiffness of the equation system is challenging for any solution strategy. The most famous and widely used solution of the last century was the algorithm of the ORIGEN [28, 29] code. The nuclei are subdivided into short- and long-living. The first ones are brought into saturation, the latter are solved for by the Taylor series expansion. This numerical scheme is fast, but has shown shortcomings in accuracy.

As already said, in MURE, two methods are actually available to solve Bateman’s equations : a 4th order Runge-Kutta method and the Chebyshev Rational Approximation Method (CRAM) which is now the default. These equations are integrated between 2 MC steps. The EvolutionControl class allows us to interact with the evolution via some methods (just after a MC run, during the evolution, just before the new MC run, see EvolutionControl class). Interactions could be done at each sub-step via EvolutionControl::ControlBeforeEachSubStep. This implies to provide the number of sub-steps in between 2 MC runs which is done via EvolutionSolver::SetNSubSteps (default is 10) or via EvolutionSolver::SetSubDT (in this case, a δt is provided, and the number of sub-steps is calculated accordingly).

7.2.3.1 4th order Runge-Kutta

The fourth-order Runge-Kutta method is a very well-known and standard method to integrate ODE ; it is generally a good compromise between speed and accuracy compare to higher (or lower) order. The method implement in MURE is adapted from Numerical Recipes[30] and is coupled to an adaptive step size. To use it:

gMURE->GetEvolutionSolver()->SetDepletionSolver("RungeKutta"); 

7.2.3.2 Chebyshev Rational Approximation Method (CRAM)

A quite elegant and accurate solution called CRAM was recently developed by M. Pusa[31, 32, 33] and has become kind of a standard solution to Bateman equations in nuclear systems. The method is very effective for matrix systems where the eigenvalues are distributed near the negative real axis. Several formulations for CRAM exist. Because of better numerical stability the method of incomplete partial fractions (IPF) is adapted with orders from N=4 up to N=48 and the CRAM coefficients are taken from the publication[33]. The IPF algorithm foresees a series of N/2 matrix inversions, which are of type sparse matrix because of the nature of the burn-up (transmutation) matrix. The matrix inversion can be done by a Gauss-Seidel iteration or – in this case – by the solver SparseLU of the Eigen library[20].

Dependent on the time range of the half-lives the recommended order lays between 16 and 48. If not requested a specific order the code chooses the suitable order for a given matrix system and time-step. CRAM has become the default solver of MURE[34]. Its options may be addressed by following methods:

gMURE->GetEvolutionSolver()->SetDepletionSolver("CRAM"); // default
// setting explicitly the order deactivates the search for the lowest suitable order
gMURE->GetEvolutionSolver()->SetCRAMOrder(48); 
gMURE->GetEvolutionSolver()->SetCRAMMethod("IPF"); // only IPF available (default)

In MURE v2 the sparse solver was based on the Eigen::SparseLU solver which shows deteriorated performance the larger the time interval becomes (in the order of 30 days). The default solver now is the Eigen internal BiCGSTAB solver. However, one may choose an UMFPACK [35] accelerated external solver of Eigen from SuiteSparse1 which could be even more powerful but it is very system dependent due to interaction of BLAS/LAPACK with pThread or OpenMP2.

7.3 Different way of evolution

We suppose that the evolution is done at a constant power. Evolution can be performed in different ways that can influence the result.

  1. Evolution with a “constant” reaction rate: this is the default. Between 2 MC steps, reaction rates are keep “constant” ; a new reaction rate evaluation is performed at the next MC run for the next step. In fact, they are not really constant: due to renormalization of the flux to keep constant power, reaction rates used in Bateman equations are in fact, ασϕ where α is the flux renormalization coefficient. Thus, periodically within a MC step, these reaction rates are readjusted (this explains the use of the term “constant” in quotes). An example of such evolution can be found in MURE/examples/Depletion/EvolvingSphere.cxx or MURE/examples/Depletion/EvolvingSphere_serpent.cxx

  2. Cross-section evolution: this is no longer the default. Previous MC runs are used to linearly extrapolate the reaction rates for the next step.

  3. Window average evolution: a mean (σϕ) is calculated from the previous MC runs and the actual one ; this mean value is used to solve Bateman’s equations. The window to compute average is based on EvolutionSolver::GetFitRangeNumber() previous runs (default is 4).

  4. Predictor-corrector methods. TO BE NOTICED: Because no real advantage of this methods have been found, this way of evolution is not more maintained. So new MURE (from MURE v2 alias SMURE) development have probably broken this way of evolution which is discouraged. There are 3 Predictor-corrector (PC) methods implemented. This is actually new and not fully tested (despite the fact that it currently works). Suppose that the MC run at time TN has been done (leading to reaction rates (σϕ)N). The evolution has to be performed up to TN+1.

    1. PC at middle point: the predictor run is performed with “constant” reaction rates (σϕ)N up to TC = T +T -N--2N+1 ; then a corrector MC run is performed and leads to (σϕ)C. Starting with the composition of TN, the evolution is done again from TN to TN+1 with the “constant” reaction rates (σϕ)C evaluated at the middle point.

    2. PC at end-point: the predictor run is performed with “constant” reaction rates (σϕ)N up to TC = TN+1 leading a composition CN+1P; then a corrector MC run is performed and leads to (σϕ)N+1. Starting with the composition of TN, the evolution is done again from TN to TN+1 with the “constant” reaction rates (σϕ)N+1 leading to a composition CN+1C. Then, the composition used for the next predictor run (for TN+1 to TN+2) is CN+1 = P C CN+1+CN+1- 2.

    3. Kind of PC at end-point: the predictor run is done like in b) but for the corrector run, we used a linear interpolation of reaction rates from (σϕ)N to (σϕ)N+1. The next composition for TN+1 to TN+2 is CN+1 = CN+1C.

    First tests seem to give an advantage to method b) for the quality of its results (compositions, dispersion of composition for M identical run with different random seeds, ...).

  5. Multi-Group evolution[11]. If this method is used, not that reaction rate tallies are not built in the MC input file ; instead, a thin energy binning flux will be built for each evolving cell. Then, after the MC run, reaction rates are calculated as ∫ ϕ(E)σr,i(E)dE where the sum is done over each energy group, σr,i is the pointwise cross-section of reaction r for nucleus i read in the ACE file used for MC run (if ASCII version is used, please note that this will take longer...). The main advantage of this method is a considerable CPU time saving (more than 30 times faster3 than the “standard” evolution). The accuracy of the result are quite good (~1% for thermal systems and 5% for fast ones4). To use this method, one has to define “gMURE->GetEvolutionSolver()->UseMultiGroupTallies();” at the beginning of the input file (see MURE/examples/Depletion/EvolvingSphere2.cxx example). The default energy binning for flux is (72995 groups)
    Energy range 10−4eV - 1eV 1eV - 10eV 10eV - 10keV 10keV - 0.1MeV 0.1MeV - 20MeV
    Number of bins/decade 2000 2000 5000 2000 1000

    The main problem of multi-group evolution is the 238U. Its capture cross section has been chosen to “optimize” the group structure (number of bins/energy interval). In order to take the advantage of the multi-group CPU time saving, it has been implemented the possibility to use multi-group runs for all nuclei but 238U : in that case, the “standard” method (real MC tally) is used for 238U, and multi-group tallies for the other nuclei. This is done by setting the “StdTallyFor238U” to “true” in the “EvolutionSolver::UseMultiGroupTallies()” method.

7.3.1 Cross-section Evolution

Cross-sections are evolving with time because of the flux shape evolution. Between 2 MCNP steps, cross-sections are either constant, either evolving. Taking the cross-section evolution during the “evolving step” should enable a better precision and increase evolution step duration, reducing the CPU time. It is based on linear extrapolation of previous MC runs.

!!!! After a study of error propagation, it appears that fitting σϕ leads to a bigger dispersion of cross-section evolution in the result (with a constant power). Thus this treatment is no more used by default. To use it you must specify it by a MURE::FitSigmaPhi(true).

7.3.1.1 How does it work?

The treatment is applied for each cross-section of each nucleus in each cell. Thus, for the sake of simplicity, we will only consider a reaction cross-section σi of a nucleus in a cell.

Since for evolution, reaction rates are used, we should rather consider the reaction rate σiϕ and not the cross-section σi. Suppose that we have performed an evolution up to step k (i.e., kth MC run is just finished). To know the reaction rate evolution up to step k + 1 (i.e., just before a new MC run), a linear fit (least square method from Numerical Recipes) is applied to EvolutionSolver::GetFitRangeNumber() previous MC reaction rates5 (that is to say, MC runs k − 3, k − 2, k − 1, and k). Then, for the evolution over step k → k + 1, the reaction rate is extrapolated with (σiϕ)(t) = at + b, where a and b are the results of the linear regression. Note that for the first EvolutionSolver::GetFitRangeNumber() MC runs, reaction rates are kept constant.

As already mentioned, integration of Bateman equations is in fact done in N sub-steps by either the CRAM method or by a 4th order Runge-Kutta steps which are automatically divided in other shorter time steps (adaptive step RK method) ; the reaction rate evolution is only performed at each sub-step level and it is kept constant within the current sub-step. This constant value is evaluated in the middle of the sub- step to be more accurate (see fig. 7.2).

pict

Figure 7.2: Illustration of a reaction rate evolution. MCNP results are in black dots. The reaction rate used during the evolution from sub-step k to step k + 1 is in green whereas the red squares indicate the time where the extrapolation is performed.

In order to better consider the rapid variation of a reaction rate slope, we have slightly improved the method as follows (here Nfit =EvoltionSolver:: GetFitRangeNumber()) :

7.4 Reactivity Control

NOTE: only few tests have been performed and therefore many bugs may be present in the Evolution Control.

As already mentioned, the control of the whole evolution is done via the EvolutionControl class. This class already includes some useful methods, but the best thing to do in order to implement your own control is to write a derived class. In a first example, the reactivity control uses the original EvolutionControl class ; then the way of implementing your own control will be given. The EvolutionControl class inheritance tree is shown on Fig. .7.3

pict

Figure 7.3: EvolutionControl class and siblings.

7.4.1 A simple example using standard EvolutionControl

The example file poison.cxx (MURE/example/poison.cxx) corresponds to a very simple reactor (cylinder of borated water with a hexahedral fuel pin lattice). The reactivity control is done by adjusting the concentration of a poison (the boron) in the water. We assume that a large reserve of boron is available and thus, when necessary, we can adjust the boron concentration at a desired value. The method principle is to calculate an estimation of keff using the reaction rate extrapolations:

-----------------------νNf-σfϕ------------------------ keff = Naσabsw∕opoisonϕ + Npσpoisonϕ− N σn,2nϕ − 2Nσn,3nϕ+ escape

(n,xn), “escape” and σpoison are supposed to be constant between 2 MC steps ; the poison concentration Np is thus determined using the desired keff in the formula.

7.4.1.1 Poison material declaration

    ControlMaterial *BoratedWtr=new Material (0.7187023, 673);
    BoratedWtr->SetDensity(0.719); 
    BoratedWtr->AddTheControlNucleus(5, 10,0.00012);
    BoratedWtr->AddTheControlNucleus(5,11,0.00048);
    BoratedWtr->AddNucleus(8,16,0.333);
    BoratedWtr->AddNucleus(1,1,0.667, "H2O");

A ControlMaterial is a Material which is used to control the evolution, and in this case, a poison to control the reactivity. The 3th an 4th lines informs us that boron 10 (and 11) are the nuclei used to control the reactivity. See the “poison.cxx” example.

This material must not be an evolving material, because we’ll consider that we can adjust boron concentration every time .

7.4.1.2 Escape calculation

In the keff formula, neutron leakage is needed ; in order to calculate these escapes, one has to define a F1 tally (current through surface) through the most outer Shape of the geometry ; in the example this leads to:

gMURE->SetOutermostShape(Vessel);

because “Vessel” is the outer cylinder of this “reactor”.

7.4.1.3 Evolution definition

You have to declare the Evolution Control object:

PoisonEvolutionControl *ECon=new PoisonEvolutionControl();
ECon->AddReactions(BoratedWtr);
gMURE->SetEvolutionControl(ECon);

The 2nd line gives the specific EvolutionControl command to add the required reactions useful for the control (poison absorption, ...)

7.4.2 Using your own EvolutionControl

The aim of this section is to help user implement its own EvolutionControl using inheritance mechanism. Then you will be able to write your own specific methods that will be called for the evolution control. Supposed we wish to create a MyEvolutionControl class: the hxx and cxx files are in MURE/example/.

NOTE: Don’t forget to redefine the EvolutionControl* Clone() as:

EvolutionControl* Clone(){return new MyEvolutionControl(*this);}

Each method and variable of EvolutionControl are known in MyEvolutionControl ; but because of the virtuality of some methods, the ControlAfterEachMCRun(), ControlBeforeEachSubStep() and ControlAfterEachEvolutionStep() methods of MyEvolutionControl will be called instead of the ones of EvolutionControl. This means that these original methods of EvolutionControl are no more called, and thus you have to read them carefully to implement what is needed or necessary for your evolution control.

Now to use this new class in Myprog.cxx:

#include <iostream> 
#include <cmath> 
using namespace std;
#include "MureHeaders.hxx" 
#include "MyEvolutionControl.hxx" 
...
int main(int argc, char** argv) 
{
    ...
    MyEvolutionControl *NewEC=new MyEvolutionControl();
    gMURE->SetEvolutionControl(NewEC); 
}

and to compile one should do

g++ -c MyEvolutionControl.cxx $MURECXXFLAG 
g++ -o Myprog Myprog.cxx MyEvolutionControl.o $MURECXXFLAG $MURELIBFLAG

where the variable MURECXXFLAG and MURELIBFLAG have been set to (bash example):

export MURE_DIR="/path_2_MURE_install_dir/"
export MURECXXFLAG="-g -fPIC -std=c++14 -I$MURE_DIR/source/include -I$MURE_DIR/source/external -I$MURE_DIR/source/thirdparty/eigen"
export MURELIBFLAG="-L$MURE_DIR/lib -lMURE -lvalerr -lmctal"

7.5 Evolution of a MCNP user defined geometry

This section is devoted to users who have already defined a complex MCNP file and would rather not use the MURE geometry capability when making evolution. Even if this task requires some work, the procedure is simple. User has to redefine in MURE style, ALL materials of its geometry and all EVOLVING cells (but not the shapes associated to the cells). Of course, user will need all files required by a standard evolution (BaseSummary.dat, AvailableReactionChart.dat, ...)

The principle is the following: a line by line of a MCNP user input file are copied, but

For material declaration, you must use the special Material constructor that takes as a single argument the material number ; this number is the one you have written in the MCNP input file. Don’t forget to declare ALL materials of the geometry (evolving or non evolving). For non evolving materials, you must ask MURE to include them in the MCNP output file by using Material::AddToGlobalNucleiVector() after the last input nucleus.

For cell declaration, you must use the special Cell constructor ; it takes the cell number (the one you have written in the MCNP input file), an evolving Material pointer and the cell volume (in m3). ONLY evolving cells should be declared.

The Original MCNP input name is defined by MURE::SetUserGeometryInputFile() method. You can also defined all desired cards of the general data block (source, physical card, ...). The rest of the evolution is completely analog to what is explained in the above sections. Don’t forget to provide a mean to calculate the tally normalization factor (you may either give a predefined power or use the MURE::SetTallyNormalizationFactor()).

An example is shown in UserGeo.cxx with the UserMCNPGeo MCNP input file.

7.6 More complex evolution conditions

Sometimes, the evolution is performed with power variations, some parameters such as boron concentration, temperature, ... are changing. This can be achieve through the EvolutionWrapper class (see EvolutionWrapper class) . Such conditions are illustrated on the following example:

EvolutionWrapper *EC=new EvolutionWrapper(); // create EvolutionWrapper
EC->AddPhase(100,3,1); // create a new 100 days phase with three MCNP steps, (1=logarithmic discretization)
EC->SetPowerLinear(0,4e4); // raise power from 0 to 40 kW
EC->SetMaterialBoronLinear(Cell_of_BoratedWater->GetMaterial(),1000e-6,500e-6); //boron depletion from 1000 to 500 ppm
EC->AddPhase(100,1,0); // a new single-step 100 days cooling phase
EC->SetPowerCooling();
EC->AddPhase(100,3,0); // a new 100 days phase with three equal(0=linear discretization) MCNP steps
EC->SetPowerConstant(4e4); // power is kept constant at 40kW in that phase
EC->Evolve(); // run the evolution (DO NOT CALL the MURE:Evolution method ; it is done here)

7.7 Equilibrium of Xe-135

In some circumstances, flux oscillations in weak coupled systems creates numerical (unphysical) Xe oscillations. These Xe oscillations enhanced the flux oscillations causing numerical instabilities. To reduce these instabilities, a special treatment has been implemented. To use this treatment, one has to call the EvolutionSolver::SetXe135Equilibrium(double teq) method, where teq is the time from which the treatment is really applied (default time value is 3 × T1∕2135 Xe ∼ 27h).

Here is the description of what is done:

Chapter 8
Looking at results from the evolution

8.1 The MURE data files

At present, results from the evolution are stored in the MC Run directory. There are 4 types of MURE files generated in addition to the MC files for each evolution step:

When using reprocessing feature DATA_bxxx or/and BDATA_bxxx files are also written just at the end of CRAM/RK sub-steps. When using PoisonEvolutionControl, additional files are also written such as POISON_PROPS which contains the poison proportion as a function of time.

8.2 Reading the data files with Tcl/Tk graphical interface scripts

This graphical interface is now deprecated. Except for ExamTree.tcl

To run these scripts from anywhere it is a good idea to add the TclScripts directory to your path. You must also set the DATADIR environment variable.

8.3 Reading the data files with ROOT graphical interface

This is a C++ GUI base on ROOT. Thus, before using this GUI, you must install a recent ROOT version (5.0 or any later version). This is easy because binary files can be downloaded from the ROOT home page (http://root.cern.ch). We have chosen ROOT because:

The GUI sources are in MURE/source/gui and the install version is in MURE/mure3/bin ; the best way to build this GUI, is to use the “build.sh --enable-gui” script and to set environment variables in your .basrc accordingly to what the build.sh script proposed. If Lapack package is available, the radiotoxicity analysis will be possible. LAPACK is a Linear Algebra PACKage used to solve the set of coupled differential equations for decays by a matrix diagonalization procedure. If this package is not installed, you have to install it (see § 1.4.4) to run the radiotoxicity module (but all other MureGui capabilities can run without Lapack). The MureGui executable has the following syntaxe:

MureGui [-help -h -onlytox -start num -font factor -color] DirName1 [DirName2 ...]

where

A module in the transport code DRAGON has been developed to copy the results in DATA ASCII file format.

The interface has a Main windows composed of 8 tabs (Inventories, Cross-sections, Fluxes, Reaction Rates, K_effective, Breeding Ratio, Radiotoxicity and Spatial Variables, see Fig. 8.1 ) and can have up to 8 plot windows (one for each tab).

Whereas MureGui is graphical interface, the “Terminal” window used to run this program is often used for some important warnings/information/output: don’t forget to look at it!

8.3.1 The 8 tabs of the main window

∑ Ntot = Ni i

pict
Figure 8.1: MureGui main window. Inventory Tab.
∑ -- -∑iϕiVi ϕ = iVi

∫ t keff(t) = 1 keff(t′)dt′ t 0

8.3.2 Frame “Cell Spatial Var” in the Main window

Users can define Spatial variables for cells (see the Cell class). The aim of these variables is only to extract results in a simpler way (this means that you can add spatial variables after the evolution is done and redo the evolution keeping MCNP files without any problems). If such variables have been defined you can plot quantities according to spatial variables.

Suppose we have done the evolution of 3 rings of a cylindrical core divided in 3 parts along the cylinder axis. If, for each of these 9 evolving cells, 2 spatial variables named “R” and “Z” with value corresponding to the typical cell locations have been defined, then the combo box of the “Cell Spatial Var” will allow users to plot for example 235U quantity for all cells with a particular radius Ri (i.e. whatever other Z value, see Fig 8.7) or to plot this quantity only for a given Ri and Zi (see Fig 8.8)

pict

Figure 8.7: MureGui - Spatial Variable combo box selected for all cells of R = 1.

pict

Figure 8.8: MureGui - Spatial Variable combo box selected for the cell located at R = 1 and Z = 0.

8.3.3 Frame “Time Interpolation” in the main window

A linear interpolation is performed on the plotted data for the selected Tab at the corresponding time (and not the selected data of this Tab). Results are returned in the terminal. Note that no time unit is specified, thus it must be coherent with the plotted data.

pict

Figure 8.9: MureGui - Interpolation and Error Bars.

8.3.4 Frame “Error Bars” in the main window

When statistical errors are available, error bars may be added to the corresponding graphs. Note that MCNP provides errors on fluxes, reaction rates and keff values. Thus, the Error Bars Button has no effect on the Inventory and Conversion Rate tabs. In [B]DATA_xxx files, are written keff, fluxes (with their MCNP absolute error), compositions of each nuclei (in atoms), and cross-sections with their error. This error is not really correct because it has been calculated as Δ-σ σ = ∘ (-----)----(--)- Δ(σϕ)-2 + Δϕ-2 (σϕ) ϕ . This is correct only if ϕ and (σϕ) are not correlated...which is obviously not the case. Thus only error for keff and fluxes are really correct.

For the average flux ϕ on several individual cells (with an individual flux of ϕi), we have:

Δϕ = -1- Vtot∘ ----------- ∑ (V Δ ϕ)2 i i i

If the sum of cross-sections is computed, the associated errors is given by:

Δσsum = -1-- Ntot -------------- ∘ ∑ (N Δ σ-)2 j j,tot j

For the average reaction rate Nj,totσϕ of one specific reaction j (linked to the corresponding nucleus) on several individual cells i, we have:

Δ(Nj,totσjϕ) = ∘ --------------- ∑i Δ(Nj,iσj,iϕi)2 = ∘ -------------[(-----)2--(---)2]- ∑i (Nj,iσj,iϕi)2 Δσσj,i + Δϕϕi j,i i

For the error associated with the power computation, we can used similar formula as for reaction rates, i.e.:

Δ(∑ jPj) = ∘ -------------------[----------------] ∑ ∑ 2 (Δ-σj,i)2 (Δϕi)2 j i(QjNj,iσj,iϕi) σj,i + ϕi

8.3.5 The Plot/Save/Quit buttons in the main window

8.3.6 More about Radiotoxicity Tab

8.3.6.1 Out flux radiotoxicity evolution

It is possible to evaluate the time evolution of radiotoxicity (or activity, ...) of an user input composition without making any MURE/MCNP evolution. This is of course an out flux evolution (only decays are taken into account). User should create a file (let say “DataFile”) in a given directory (let say “ToxDir”). This file must contain

V -1
UNIT
n
Z1 A1 I1 M1 P1
Z2 A2 I2 M2 P2
...
Zn An In Mn Pn

where UNIT is a string among “At”, “Mol”, “Wgt(g)”, “Wgt(kg)”, “Wgt(t)” giving the unit of each “proportion” Pi ; “n” is the number of different isotopes (Zi, Ai, Ii) and “Mi” is the molar mass of the isotope i. For example, to see the evolution of 2.975 kg of 238U the U8_data file is

V -1
Wgt(kg)
1
92 238 0 238.05 2.975

Then launch MureGui as :

MureGui -onlytox ToxDir/U8_data

and go in the “Radiotoxicity” tab ; this gives a total radiotoxicity of about 1mCi for 238U at t = 0 ; then it reaches the secular equilibrium after 2 million years to 14 × 1mCi (13 radioactive daughters of 238U+238U).

Remarks:

8.3.6.2 Gamma, Beta, Alpha, and Neutron spectra of an evolving material

User can plot spectra of an evolving material. This feature is “hidden” in the Radiotoxicity Tab : press the “CHARS logo” or depending on your installation, the small down arrow, between “Irradiation” and “Cell Spatial Var” to make it appear (fig. 8.10) .

pict

Figure 8.10: MureGui Radiotoxixity Tab.

First of all, spectra has to be built pressing the “Build Spectra” button. It will create a γ spectra for the “Gamma Spectrum” tab, β, α and neutron spectra in the corresponding tabs for each nucleus. This step could take a while especially for the Gammas. In the “Gamma Spectrum” tab (see fig. 8.10), gammas produced by decay modes selected in the “Decay Mode” from are plotted.

“Plot Spectrum” frame: This frame is used to plot either 3D spectra (activity versus energy and time) for γ, β, α and neutrons or “time slice cuts” using the “2D slice” slider of these spectra. For 3D spectra, only one nucleus (or one group among “total”, “FP” or “MA”) can be plotted.

“Find emitter” frame: This frame (only for gammas and alphas) allows user to find the emitter of rays of energy E ± ΔE enters in combo boxes (E and ΔE are both in keV, see fig. 8.11 ).

pict

Figure 8.11: Find gamma emitters at 185keV ± 1keV . They are sorted by decreasing intensity (gamma intensity×atoms) and the decay mode is indicated (alpha, beta, ...).

“Save Data” button: This button can save the last plotted values in either ASCII format or “MCNPSource Input” format (only working for “2D Slice” plot)

“Use binning from file” button: This button is used to load the energy binning for Spectrum from a file. Warning: the starting reading path for the file is $MURE_DIR/share/data/GuiData/.

The format of this file is the same of the one used for Spectrum class (see for example $MURE_DIR/share/data/GuiData/NEA_Binning_Gamma – energies are in eV):

NumberOfBins
LowerEnergy[bin=0] //lower bound of energy bin in eV
LowerEnergy[1]
...
LowerEnergy[NumberOfBins]
UpperEnergy[NumberOfBins+1] // upper bound of the last energy bin in eV

If a spectrum of a given cell or of sum of selected cells is wanted, check the wanted cell(s) in Flux Tab.

8.3.6.3 Radiotoxicity of a sample

To compute inventories, decay heat, radiotoxicity and alpha, beta, gamma and neutron spectra of a sample over the time one can use the option “-onlytox” of MureGui.

First define the sample Material :

Material *MySample=new Material();
MySample->SetDensity(Density);
MySample->AddNucleus(Z,A,I,Atomes);
...

and then dump it:

MySample->DumpMaterial(“SampleName”,Mass,Volume); //Mass [kg], Volume [m3]

Finally open it with MureGui :

Write in you terminal :

MureGui -onlytox SampleName

8.4 Spectrum Classes (MCNP only)

The aim of these classes is to characterize any irradiated sample in order to create elaborate sources for radio-protection studies. Spectrum class is the base class and have not to be used directly bu via its daughters (GammaSpectrum, NeutronSpectrum, …)

8.4.1 Spectrum Class

It’s a base class to define energy spectra for γ, β, α and neutrons. This class defines Spectrum objects, i.e. one array store intensities an other store energies. The size of these array are NumberOfBins+1. To declare a Spectrum proceed as follow :

Spectrum* MySpectrum=new Spectrum(EnergyMin,EnergyMax,NumberOfBins); //Energy unit: eV

Creates a energy spectrum with a constant linear binning from EnergyMin to EnergyMax (see Spectrum for more options).

For non constant binning, one can use:

Spectrum* MySpectrum=new Spectrum(EnergyLowVector);

where EnergyLowVector is the vector of the lower bounds of each energy bin (NumberOfBins elements) except the last element (the (NumberOfBins+1) of the vector which is the upper bound energy of the last bin. Other methods allow us to build spectrum binnings from arrays or from a file.

To build a spectrum for a particular radiation one uses the GammaSpectrum, BetaSpectrum, AlphaSpectrum or NeutronSpectrum classes :

XSpectrum* MySpectrum=new XSpectrum(EnergyLowVector); //where X can be Alpha, Beta, Neutron or Gamma.
MySpectrum->Fill(MyMaterial,Volume);

This last method compute the alpha, beta, gamma or neutron spectrum of MyMaterial of volume Volume.

A spectrum can be added to an other with Spectrum::AddSpectrum() but only if the two spectra have the same energy binning:

MySpectrum->AddSpectrum(SpectrumToAdd);

If you add a XSpectrum with a YSpectrum, the particle type, that can be used for an MCNP source, is defined by the one of the object calling the method (MyGammaSpectrum->AddSpectrum(MyNeutronSpectrum) will have for particle type γ).

A Spectrum can be dump in a file with (if WithEnergies is true, the energy bounds of each bin is given, else (default), the middle energy of the bin is given) :

MySpectrum->Dump(filename,WithEnergies);

All the following classes inherit from the Spectrum class: AlphaSpectrum, BetaSpectrum, GammaSpectrum and NeutronSpectrum.

8.4.2 GammaSpectrum class

In the method GammaSpectrum::ReadENSDF(Z, A, I, Atoms, DecayMode) where Atoms is the number of atoms of the ZAI and DecayMode is the decay mode leading to the emission of a gamma ray (it is any combination of “A” (α decay), “B” (β− decay), “E” (electronic capture) and “I” (Isomeric Transition).

GammaSpectrum::ReadENSDF() searches in ENSDF files the gamma transition of theZAX daughter(s) from these values : decay constant of ZAX (λ ), Branching Ratio of the decay (BR), Relative Intensity (RI), and the Normalization Factor (NR). It calculates the activity A for one transition of Energy E and fills the Spectrum:

A = λ× NA × RI × N R × BR ∕100[bq] ZX

where RI ×NR is the number of photon per 100 decays of the parent for this decay branch. RI ×NR×BR is the number of photon per 100 decays of the parent. In order to compute the gamma spectrum of a sample, have a look in MURE/example/SpectrumExample.cxx.

8.4.3 AlphaSpectrum class

Here, the method AlphaSpectrum::ReadENSDF(Z, A, I, Atoms) is almost identical to the GammaSpectrum::ReadENSDF() one, but it looks for alpha ray energies and intensities instead of gammas.

8.4.4 BetaSpectrum class

The method BetaSpectrum::ReadENSDF(Z, A, I, Atoms) finds in ENSDF files decay constant, branching ratio, normalization factor, intensity, and also the the endpoint energy (or data needed to calculate it if it’s not present, i.e. energy of the decaying level, Q-value (ground state to ground state) of the decay and energy of the level of the daughter)

EndP ointEnergy = Qvalue+ EnergyOf TheDecayingLevel − EnergyOf TheLevelOfT heDaughter

As soon as EndPointEnergy is known, the distribution f(E) is calculated using a simplified Fermi theory where all the transition are allowed, then for each bin, the Intensity is multiplied by this normalized distribution.

A(E) = Atoms × λ × RI × NR × BR ×f (E )

where

∘ ------- f(E) = Norm × ----2Πη------× W × W 2 − 1 ×(EndpointEnergy − E)2 1− exp(− 2Πη)

with η = 2 Zdau4gΠh𝜀0teℏrV×ee-, V e = c ×∘ ---------- 1− (1+1W)2 and W = 1 + meEc2 .

8.4.5 NeutronSpectrum class

This class contains methods which are required to calculate neutron energy spectra. Neutrons taken into account are neutrons from spontaneous fission (NeutronSpectrum::ReadSFData method), neutrons from (α, n) reactions (NeutronSpectrum::AddAlphaNSpectra method), and neutrons from β−N (NeutronSpectrum::ReadENSDF method).

8.4.5.1 Neutron from spontaneous fission

The neutron spectra from spontaneous fission is computed using a normalized Watt distribution W(E)

W(E) = Norm.exp(− E) a.sinh( ) √b.E-where a and b are coefficient depending on the nucleus.

Finally the neutron activity of the nucleus ZAX by spontaneous fission at energy E is :

A(E) = λ × NA X × ⟨ν⟩× BR × W (E ) Z

where BR is the branching ratio of the spontaneous fission and ⟨ν⟩ the average neutron per fission. The parameters ⟨ν⟩, BR, a and b are available in literature and those used in MURE are located in /Path_2_Mure/MURE/data/SFNeutronSpectraData.dat [36].

8.4.5.2 Neutron from (α,n) reactions

Neutrons can be produced by (α,n) reactions on light nuclei during slowing down of alpha particle in the sample. Here the reactions taken into account are only 17O(α,n)20Ne and 18O(α,n)21Ne. To compute neutron spectra from these reactions, one needs alpha spectra, stopping power in the media, cross sections leading to different level of ANe and oxygen density (for each isotopes). Default values are natural oxygen density in UOx fuel (ρ = 10.4g.cm−3).

The stopping power was calculated for a UO2 media using SRIM[37], (α,n) cross section are from JENDL(JENDL (α,n) Reaction Data File 2005 (JENDL/AN-2005)). The probability pi of an alpha of energy Eαhaving a (α,n) reactions on target i while alpha slowing down on dx is pi = Niσidx, where σi is the (α,n) cross section on target i. The probability pi can be expressed as a function of the variation of energy of the alpha particle on dx : pi = NσddExdE. By integrating the energy of the alpha particle from it’s initial energy to 0:

∫ Eα Pi(Eα) = Ni σi(EdE-)dE 0 − dx

Assuming that neutrons are emitted isotropically in the center of mass, neutrons are uniformly distributed in energy in the laboratory system between a maximum energy and a minimum energy defined by the conservation of momentum and energy (see for example [38] for more information).

8.4.5.3 Neutron from (β−n) decays

8.4.6 Define MCNP Source with Spectrum object

A computed Spectrum can be used to define a MCNP source (see also section 5.2.3, and examples SpectrumExample.cxx, and TubeSource2.cxx in MURE/example).

MCNPSource *MySource=new MCNPSource(NumberOfParticle);
MySpectrum->Fill(MyMaterial,Volume);
MySource->UseThisEnergyDistribution(MySpectrum,DistributionNumber);
gMURE->SetSource(MySource);

where Volume is the Volume of your sample and DistributionNumber is the MCNP source distribution number (default 800).

Chapter 9
Thermal-hydraulics/neutronics coupling

9.1 Preliminary Remarks

In MURE 2/SMURE, some classes have been removed (ReactorChannel) or renamed (ReactorMesh) and completely rewritten. There are no more backward compatibility of the renamed classes with older ones. Read carefully the following guide and the provided examples.

Coupling neutronics and Thermohydraulics (like the examples given in MURE/examples) is relevant only if temperature dependent cross-sections for MC transport codes have been provided for each isotope.

9.2 General Overview

Simulations of nuclear power plants can need a discretization of the core in 3 dimensions. The ReactorAssembly classes provide an easy way to build a full core (with radial and axial discretization) for hexagonal or cubuoid assemblies. These classes allow to define automatically radial zones and axial levels per assembly from a general pins (built via the PinCell class) scheme. In each radial zone and axial level, all pins of the same type (e.g. Fuel pins, Guide Tube pins, ...) will have the same composition (all fuel pins have the same clad, coolant and fuel composition in the radial zone, idem for Guide Tube, ...). Zoning is defined via virtual circle (center at in the middle of the assembly) or via given explicit positions. For the former case, a pin belong to a zone if its most outer layer (e.g. clad radius) is fully inside the given circle. These assemblies can then be used to fill a core.

COBRA_EN class has been created for the coupling analysis with a qualified thermal-hydraulics sub-channel code : COBRA-EN[14, 15] (Coolant Boiling in Rod Arrays). The coupling is “transparent” for users: the ReactorAssembly is passed to COBRA_EN class. At the present time, COBRA_EN class only support cuboid ReactorAssembly with the following restrictions:

BATH solves the heat equation in case of a permanent regime and a cylindrical fuel and cladding shapes. The following estimates should be taken into account:

9.3 Description of methods

9.3.1 Introduction

A nuclear reactor is mainly controlled by two disciplines: neutronics and thermal-hydraulics. In fact, reactor core temperatures depend on the heat sources, and thus on the distribution of power which evolves in the course of time (calculated by neutronics codes). Yet, these codes require cross sections which depend on the temperatures, and thus on the nature of flow that is solved by thermal-hydraulics codes.

The input file for a coupled analysis is almost the same than a pure neutronics MURE calculation. Of course, some thermal-hydraulic characteristics must be added (like coolant temperature entrance, pressure,...) but a non thermal-hydraulic physicist can quickly perform a simulation with a thermal-hydraulics coupling. The generation of the COBRA input file happens automatically and without any user intervention as well as the coupling procedure between MURE and COBRA.

The neutronics code calculates the power distribution in the fuel pins, which is transferred to the thermal-hydraulics code. Then, this latter determines new properties in each pin and sub-channel (temperatures, densities,...). Finally, these data are used as new inputs for the next neutronics calculation. Iterations can be repeated until a converged state is obtained or during a step time for an transient analysis.

Meshing grids are different between neutronics and thermal-hydraulics (c.f. figure 9.1).

pict

Figure 9.1: Meshing grids.

In fact two different physics disciplines are dealt with, but their solving methods have nothing in common. This is the reason why two different meshes are defined. Consequently classes manage the correspondence and the symmetrization of these meshes. Updates are done using mean values by zones (c.f. figure 9.2).

pict

Figure 9.2: 4 Radial zones and 10 axial levels.

These zones are chosen by the user and are dependent on the precision required. In each radial zone and axial level, neutron flux, reaction rates, properties of materials or coolant and fuel burn-up are treated separately.

Any increase in the discretization must be carefully considered because it may use more computing resources for a marginal accuracy gain, or even, in certain circumstances, cause a loss of precision.

9.3.2 ReactorAssembly class

The ReactorAssembly class is a complete rewriting of ReactorMesh class. This class can be used in a purely neutronics way or for thermohydraulics/Thermics-neutronics coupling. As already emphasis, when COBRA_EN is involved, some restrictions are mandatory (see § 9.2).

Utilization for pure neutronics calculation

Examples are provides in RA_1ass.cxx (1 hexagonal assembly), RA_1ass_serpent.cxx (1 cuboid assembly) and RA_lattice.cxx (a cylindrical core filled by cuboid assembly).

IMPORTANT : The order of calling method in the ReactorAssembly class is not arbitrary! This is because some methods defined and/or calculate some quantities depending on previous ones. This normally has been protected (by a MURE error). But to avoid problems here are some rules (methods not mentioned there after can be called without specific order):

  1. Define a ReactorAssembly object via ReactorAssembly::ReactorAssembly(double, int, int, string)

  2. Define eventually plenum thickness and materials if they exist via ReactorAssembly::SetPlenum(double, Material*, Material*) and assembly duct (assembly channel) if needed via ReactorAssembly::SetAssemblyDuct(double, Material*)

  3. Assign the assembly Shape to the object via ReactorAssembly::SetAssemblyShape(Shape_ptr )

  4. Add all standard fuel pins that fill the whole assembly (without specify position) via
    ReactorAssembly::AddFuelPin(PinCell*)

  5. Add specific pins at the given position via ReactorAssembly::AddFuelPin(PinCell*, int, int) or via ReactorAssembly::AddGuideTube(PinCell*, int, int) or via ReactorAssembly::AddControlRod(PinCell*, int, int)

  6. ...

  7. and last method to be called1 must be ReactorAssembly::BuildAssemblyGeometry()

As an example, user first defines Materials and PinCell (Fuel, Guide Tubes, Control Rods, ...) he want to use. Then, he defines the shape of the assembly (either a cuboid or a hexagon):

Shape_ptr AssemblyShape(new Brick(S/2., S/2., H/2.)); //cuboid assembly S x S x H

or

Shape_ptr AssemblyShape1(new Hexagon(H/2., S)); //Hexagon assembly with edges of size S and height H.

This shape determines the type of the ReactorAssembly (Cuboid or Hexagonal) ; if, for some reasons, the assembly shape is more complex (e.g. a cuboid assembly with a hole inside that is not filled by any pins, user should provide this complex shape as well as the real cuboid (Brick class) or hexagonal (Hexagon class) of the outer bounding shape).

Then user defines a ReactorAssembly object by given the pitch of the pin lattice and the number of radial zone and axial levels

ReactorAssembly* RA = new ReactorAssembly(pitch,4,10); // 4 radials zones, 10 axial levels

If equal size upper and lower plenums are needed, user should call ReactorAssembly::SetPlenum(thickness, PlenumMaterial) method before any other call to a ReactorAssembly ’s method.

The shape of the assembly is passed to the ReactorAssembly object:

RA->SetAssemblyShape(AssemblyShape); // physical limits

And zoning is defined by given as many radial zones as it has been defined in the constructor (in this example, 3 radial zones) ; zone number range should be 0 for the most inner zone to N for the most outer. To be noticed: all lattice cell of an assembly (i.e. pins) must belong to a defined radial zone ; this implies that the last zone must include all the pins. If zones are defined via concentric virtual circle center in the middle of the assembly, the last zone radius must be at least equal to S∕√ - 2 for a cuboid assembly of side S and to 2 × S for an hexagonal assembly of side S.

// Automatic differentiation of zones for structurally identical elements 
//                (c.f. fig 9.2)
RA->AddCircleZoneRadius(2.5*pitch,0); 
RA->AddCircleZoneRadius(5*pitch,1);
RA->AddCircleZoneRadius(8*pitch,2);
RA->AddCircleZoneRadius(S/sqrt(2),3);

Then, add the pins at a given position in the (x,y) lattice : the position is a position index in the lattice (with the respect to the primitive cell lattice) ; thus the (0,0) correspond to the central primitive cell, (-1,-1) to the left down adjacent cell, ... If no position is provided, the the pin fill the whole lattice (the normal way of filling normal fuel pins) ; then user can override some of the pins by adding some specific ones to a given position

RA->AddFuelPin(FuelPin); // the FuelPin (PinCell *) pins fill the whole assembly.
RA->AddGuideTube(GuideTube,5,2); // add a guide tube pin at (5,2) from the assembly lattice center
RA->AddGuideTube(GuideTube,0,0); // add a guide tube pin in the center of the assembly
RA->AddGuideTube(GuideTube,-5,-2);// add a guide tube pin at (-5,-2) from the assembly lattice center
...

You can add different type of fuel pin, control rods (with ReactorAssembly::AddControlRod), ...

And then build the assembly geometry (this MUST be the last method of ReactorAssembly called):

RA->BuildAssemblyGeometry();

Clearly, the outside cell of the assembly cannot be defined automatically and it should be defined by user himself.

For COBRA-EN coupling, the coolant around guide tubes, control rods and fuel pin cells is the same as the ones for neighbours cells (same composition, temperature, density, ...) of the same radial zone.

A more complex example is presented in RA_complex.cxx (c.f. figure 9.3) : this example shows how to define radial zones using virtual circle radii as well as specifying the zone by a lattice position.

pict

Figure 9.3: A complex assembly with 7 radial zones: 3 zones for standard fuel pins, 3 zones for specific fuel pins and 3 zones for guide tubes..
9.3.2.1 Using a ReactorAssembly to fill a core lattice

After defining a ReactorAssembly object with its fuel pins, guide tubes, plenums, ..., this assembly can be used to fill a reactor core. In order to do it, before calling the ReactorAssembly::BuildAssemblyGeometry(), one has to call the ReactorAssembly::SetUniverse() method. Then, one as to used the assembly as a lattice generator and to avoid problem of chain surface, the ReactorAssembly shape is taken as an infinite one:

...
ReactorAssembly* RA = new ReactorAssembly(pitch,3,3);
...
RA->SetUniverse();
RA->BuildAssemblyGeometry();
Shape_ptr InfinitAss=RA->GetInfiniteAssemblyShape();
InfinitAss->SetUniverse();
 
InfinitAss>>Core;
 
LatticeCell* LatticeGenerator=new LatticeCell(InfinitAss);
LatticeGenerator->FillLattice(RA->GetUniverse())

In this short example, the fuel core is filled with the assembly, but a more realistic filling scheme is presented in MURE/examples/RA_lattice.cxx (see Fig. 9.4).

pict

Figure 9.4: A core filled by a ReactorAssembly: Axial (rigth) and radial cut (left) ; a zoom on the ReactorAssembly object is also shown (middle).

For evoltion example, see the documented example MURE/example/Evolution/RMAssembly.cxx (c.f. figure 9.5).

pict

Figure 9.5: Example of an assembly created with ReactorAssembly.

9.3.3 The COBRA_EN class : coupled neutronics/thermal-hydraulics calculations

This class manages the Oak Ridge National Laboratory code COBRA-EN (COolant Boiling in Rod Arrays). It is a sub-channel code that allows steady-state and transient analysis of the coolant in rod arrays. The simulation of flow is based on a three or four partial differential equations: conservation of mass, energy and momentum vector for the water liquid/vapor mixture (optionally a fourth equation can be added which tracks the vapor mass separately). The heat transfer model is featured by a full boiling curve, comprising the basic heat transfer regimes: single phase forced convection, sub-cooled nucleate boiling, saturated nucleate boiling, transition and film boiling. Heat conduction in the fuel and the cladding are calculated using the balance equation.

Some modifications on the source code have been done to optimize calculations on Intel Fortran Compilers (ifort) and to enlarge arrays (parameters IDSIZE in source files).

9.3.3.1 Coupling Thermo-hydraulics and neutronics

  1. In the MURE input file, you must provide a gMURE->GetTemperatureMap()->SetDeltaTPrecision(Delta_T); where Delta_T is the temperature precision (in K) used to choose the closest temperature for cross-sections (this number could be big, let say arround 1500K).

  2. Then, each material where temperature evolution is wanted should know it via the Material::SetTemperatureEvolution() ; this should be done for the fuel as well as for the coolant part.

  3. a ReactorAssembly is built (as described above but be sure to have created real clones of evolving temperature’s materials ; for example the coolant used in fuel pins, plenums, guide tubes,... should be a different material.

  4. Create a COBRA_EN object (let say this is a pointer called Channel) with the previous ReactorAssembly.

  5. Set to this object the boundary conditions (inlet temperature, mass flow, exiting pressure,...), see § 9.3.3.2

  6. then build and run the COBRA-EN file

  7. Optionally, to achieve temperature convergence, one can loop on successive iterations between MC/thermo-hydraulics calculation.

All theses steps are exemplified in MURE/examples/Thermal/SimpleTHAssembly.cxx.

9.3.3.2 Input/Output file

The run of successive MC code and COBRA code required input files for both codes and transfer files that allow communication and cross-check of the simulation. MC file (inpXXX) is generated as usual in MURE, in the (MC) run directory. In this directory, a “TH_yyy” (where yyy is the name of the reactor assembly), input/output and transfer file for COBRA-EN are written.

Input data on operating conditions

Some precision about fluid mechanics must be given. Otherwise, default values are taken into account.

Inlet temperature (K) :

Channel->SetInletTemperature(560);

Average inlet mass flux (kg∕s∕m2) :

Channel->SetInletMassFlux(2923.72);

Inlet boron concentration (ppm) :

Channel->SetInletBoronConcentration(1000);

System exit pressure (MPa) :

Channel->SetExitPressure(15.8);

See the documented examples MURE/examples/Thermal/SimpleTHAssembly.cxx and MURE/example/Thermal/C4Assemblies.cxx (c.f. figure 9.6) and some results :

pict
Figure 9.6: Geometry visualization.
pict pict
Figure 9.7: Complex COBRA coupled simulation : 20 axial levels and 4*7 radial regions. Right: 3D Power distribution (W). Left: Exit coolant temperature (K)

9.3.4 The BATH class (NOT YET TESTED IN SMURE-MURE v2.0)

This class manages temperature calculation for a coupled simulation neutronics/thermal-hydraulics with a 2D RZ precision. It is not as accurate as COBRA and must be explicitly coupled by the user but its advantages are :

9.3.4.1 Capabilities

The code already includes data and equations necessary for simulations. It uses:

9.3.4.2 What is solved

9.3.4.3 How to add data

If thermal conductivity of new cladding is needed, one has to modify the source of ThermalCoupling class, adding flags and switching cases in ThermalCoupling::CladConduction().

To add a new coolant, equations of dimensionless number calculation are perhaps necessary and data tables must be added in thermal_data directory. As it is done for water and sodium, files containing thermal data must be created : see Water.xls in thermal_data directory.

Adding a new fuel is easier: only 2 tables of data are necessary (density and thermal conductivity : see thermal_data directory).

Contact Nicolas Capellan (capellan at lpsc.in2p3.fr) for further explanation.

9.3.4.4 Use the code

Users should create the geometry themselves. Flags will be added to indicate positions of frames and cells to the thermal-hydraulics code :

BATH is created and launched after the MCNP run.

See the documented example MURE/examples/Thermal/BATHExample.cxx.

9.3.4.5 Output data

Outputs are created in the directory called Zi_Run_j (i is the “region” and j the iteration number of coupled field).

For each level, thermal results are built together with the temperature function of distance from the center (temperature gradient inside fuel, cladding and coolant, c.f. figure 9.8).

pict

Figure 9.8: BATH output example.

Appendix A
Basic of C++ to understand MURE

A.1 Introduction

This appendix is not a C++ tutorial ; its aim is to define roughly commonly used words inside MURE user guide and to allow users unfamiliar with C/C++ to understand examples and uses of MURE. Each main program in C++ is defined as

int main() //or int main(int argc, char **argv) to use argument passed to the program
{
   ...
}

In general, special functions, variables, ... are needed. Then, the definition of these functions known as the header files should be included:

#include <iostream> //standard functions for input/output on the screen (cout, cin, ...)
#include <cmath> //standard mathematical function (cos, sin, ...)
using namespace std; 
#include “my_header_file.hex” //a definition made by a user
 
int main()
{
...
}

A.2 Class & Object

C++ is said Object Oriented. This means that the notion of variables is extended to a more complete and self consistent form: the object. The variable types are now called classes: for example, the type float of C (real in fortran) is now a class. For this class the difference between type and class is simply semantic. But you can create your own type (your own class) that regroups many other objects (or variables in C/fortran language) of different classes (different types). For example a class Point can have 2 coordinates of the class double:

class Point
{
   public:
     double X;
     double Y;
};

The “public:” keyword means that one can use the objects of this class from the outside of the class. The 2 objects X and Y are called members or attributes of the class Point. To declare an object P of the class Point to be a point of coordinates (0,4), one has to

Point P;
P.X=0;
P.Y=4;

Moreover, a class contains methods, that is to say, specific functions that can act only on objects of this class:

class Point
{
   public:
     void Translate(double dx, double dy) {X+=dx;Y+=dy;}
     double X;
     double Y;
};

P can be translated of V (1,1):

P.Translate(1,1);

Special methods of classes are named constructors: these methods have no type and name of the class ; they are used to allocate memory and eventually initialize an object:

class Point
{
   public:
     Point(double x, double y){X=x; Y=y;}
     void Translate(double dx, double dy){X+=dx;Y+=dy;}
     double X;
     double Y;
};

Then P can be declared as

Point P(0,4);

Methods (as well as functions) can be overloaded: 2 methods may have the same name but different arguments:

class Point
{
   public:
     Point(double x, double y){X=x; Y=y;}
     Point(){X=0; Y=0;}
     void Translate(double dx, double dy){X+=dx;Y+=dy;}
     void Translate(double dV[2]){X+=dV[0];Y+=dV[1];}
     double X;
     double Y;
};

here the second constructor has no argument ; it is called the “default constructor” ; in this case it defines the origin point ; the second Point::Translate() method needs an array of dimension 2 as argument:

Point P1; //call the default constructor
Point P2(0,4); // call the “normal” constructor P2 is (0,4)
double dv[2]={1,1}; //an array of 2 doubles
P1.X=2; //set X of P1
P1.Y=2; //set Y of P1;
Point P3=P1; //P3 is (2,2)
P1.Translate(dv); //Translate P1 from (2,2)->(3,3) but P3 is still at (2,2)
P2.Translate(1,1);//Translate P2 from (0,4)->(1,5)

It is usually safe to “protect” the objects of a class from modifying them from the outside of the class without using a class method. One may succeed by using the “private:” keyword. But then one must define new methods to modify them.

class Point
{
   public:
     Point(double x, double y){X=x; Y=y;}
     Point(){X=0; Y=0;}
     void Translate(double dx, double dy){X+=dx;Y+=dy;}
     void Translate(double dV[2]){X+=dV[0];Y+=dV[1];}
     void SetX(double x){X=x;}
     void SetY(double y){Y=y;}
     double GetX(){return X;}
     double GetY(){return Y;}
  private:
     double X;
     double Y;
};

Then, one should not write

P.X=1;
cout<<”Coordinate of P(”<<P.X<<”,”<<P.Y<<”)”<<endl;

but instead, use

P.SetX(1);
cout<<”Coordinate of P(”<<P.GetX()<<”,”<<P.GetY()<<”)”<<endl;

Of course, this is not very useful for this simple example ; but for larger codes, it avoids bugs.

A.2.1 Header and implementation files

In the previous class Point, one has defined explicitly all methods inside the class. But when dealing with larger classes, libraries, please consider separating the class definition (.hxx files called headers) from the class implementation (.cxx files). Then, one may build libraries with implementation files whereas the header files could be included in any file using the library. For example one may rewrite the Point class header in the Point.hxx file

class Point
{
   public:
     Point(double x, double y);
     Point();
     void Translate(double dx, double dy);
     void Translate(double dV[2]);
     void SetX(double x){X=x;} // here the methods are very short, 
     void SetY(double y){Y=y;} // one can let then with inline definition.
     double GetX(){return X;}
     double GetY(){return Y;}
  private:
     double X;
     double Y;
};

Not the “;” after the first 4 methods: it means that they will be defined somewhere else. The Point.cxx implementation file:

#include “Point.hxx”
Point::Point(double x, double y)
{
   X=x; 
   Y=y;
}
Point::Point()
{
   X=0;
   Y=0;
}
void Point::Translate(double dx, double dy)
{
   X+=dx;
   Y+=dy;
}
void Point::Translate(double dV[2])
{
   X+=dV[0];
   Y+=dV[1];
}

Note the “Point::” ; it means that the following method belongs to the class point (“::” is the scope resolution operator).

A.3 Default arguments

In MURE, lots of parameters have default values. In general default values are assigned in methods (or function) definitions (header files) with the following syntax:

void TheMethod(double x=3);

This means that if one calls TheMethod(5), the value of the argument x=5, whereas a call to TheMethod() will take x=3 for this “default” argument. If one writes

void AnOtherMethod(double x, double y=2);

The first argument must be given and the second one (y), if it is not given, takes the value y=2. All default arguments must be defined after the explicit argument1. If one writes

void ThirdMethod(double x=1, double y=2, double z=1);

then a call to

ThirdMethod();         // means that x=1, y=2, z=1
ThirdMethod(3);        // means that x=3, y=2, z=1
ThirdMethod(3,10);     // means that x=3, y=10, z=1
ThirdMethod(3,10,100); // means that x=3, y=10, z=100

A.4 Pointers

We have defined objects of classes. One can define pointers on objects:

Point *P1=new Point; //call the default constructor
Point *P2=new Point(3,4); //call the normal constructor
P1->SetX(1); //because P1 is a pointer on a Point, one use “->” instead of the “.” 
P2->Translate(1,1);
Point *P3;
P3=P1; //P3 and P1 are pointing on the same Point (1,0)
P3->Translate(1,1); // P3 and P1 point on (1,1).
Point *P;
P=new Point[3]; //P is an array of 3 Points (it points on the first)
for(int i=0; i<3; i++)
{
   P[i].SetX(2*i+1); //P[i] is a Point (whereas P is a pointer)
   P[i].SetY(i+1);   //thus one uses a “.” and not the “->”
}
P[0].Translate(1,1); //Translate P[0] from (1,1)->(2,2)
delete P1;
delete P2;
delete [] P;

The 2 first lines allocate 2 Points (the “new” keyword). One use the “->” instead of the “.” to use a method (or a public attribute) because P1 and P2 are pointers. In the declaration of P3, no Point is allocated (no “new”) ; the “P3=P1” line means that P3 points on the same Point than P1 ; thus translating P3 will translate P1 (and vice verse). By writing “P=new Point[3]”, one allocates an array of 3 Points by calling the default constructor of Point (to declare arrays, a default constructor must exist). Each P[i] is a real Point object (not a pointer). When one use pointers, the memory is not freed by itself: one has to do it manually though a “delete”. Here one can delete P1 or P3 but not both ; in the example, P3 continues to point on the address of the P1 Point, but because the memory has been freed, this address does not correspond to anything. To delete arrays, one has to use the “[]” after the “delete” keyword.

A.5 Inheritance

One very powerful tools of classes is the inheritance mechanisms. One says that a class inherits from a mother class: all members (methods and attributes) of the mother class are known from the Daughter class. This one may have new methods and/or attributes. Following our previous example, one can define a Point3D as:

#include “Point.hxx”
class Point3D : public Point //here one says that Point3D is the daughter of Point
{
   public:
     Point3D():Point(){Z=0;}
     Point3D(double x, double y, double z):Point(x,y){Z=z;}
     void Translate(double dx, double dy, double dz){X+=dx;Y+=dy; Z+=dz;} //possible if the “private:” of Point
     void SetZ(double z){Z=z;}                                            //is replaced by “protected:”
     double GetZ(){return Z;}
   private:
     double Z;
};

and one has to replace the keyword “private:” of the class Point by the keyword “protected:” in order to allow the Daughter of Point (Point3D) to modify its attributes. But from the outside of Point and Point3D, it is impossible to modify these attributes. The 2 constructors of Point3D start to call respectively the default and the normal constructor of the mother class (:Point() and :Point(x,y)).

Point3D P1(1,2,3); // P1 is at (1,2,3)
Point3D P2; //P2 is at (0,0,0)
P2.SetX(1);
P2.SetY(1);
P2.SetZ(1); //P2 is now at (1,1,1)
P1.Translate(1,1,1); //P2 is now at (2,2,2);

A.6 Namespace

Some classes or methods are defined in name spaces : this means that one can use a class/method with a reference to that name space : for example the use of “cout” to print a flow should be done as

#include <iostream>
int main()
{
    ...
    std::cout<<”bla bla bla”<<std::endl;
    ...
}

But because this is sometimes quite heavy to write, one just as to use

#include <iostream>
using namespace std;
int main()
{
    ...
    cout<<”bla bla bla”<<endl;
    ...
}

This trick is also used in MURE to choose the kind of output one need (Serpent or MCNP) ; lots of classes have the same name (e.g. the Brick class) but, because they belong to 2 different name spaces, they are doing different things. Thus one can use either MCNP::Brick or Serpent::Brick or more easily, for a “Serpent Brick”,

#include <iostream>
using namespace std;
using namespace Serpent;
int main()
{
    ...
    Brick *b=new Brick(...);
    ...
}

Appendix B
Node tree simplification

In this appendix, more details are provided on the principle of the Node simplification to obtain more optimized MCNP input files. For example, in figure B.1 the Plane number 2 belongs to 2 Nodes and must not be destroyed after the destruction of the included node.

pict

Figure B.1: A tree corresponding to a geometrical object before its simplification. The representation tells whether a node is a union or intersection, then the node address is given, the address of the not of the node (complementary node) is also given if it exists. Nodes in red are usually not removable because of included, bounding or inside shapes (see Node class). For each leaf of the tree, the type of Shape (Plane, Sphere, ...) is given as well as its address and its surface number with a sign (like in surface in MCNP cell).

pict

Figure B.2: The tree of figure B.1 after simplification. Since one node was included in the other, it has been removed.

pict

Figure B.3: An other example of tree before simplification.

pict

Figure B.4: The tree of figure B.3 after simplification.

Appendix C
Stochastic volume calculation

WARNING: UP TO NOW SURFACE CALCULATIONS HAVE NOT YET BEEN IMPLEMENTED.

When using Tallies, cell volumes and/or surfaces are needed. In this package, cell volumes are calculated if possible ; when the program is not able to calculate cell volumes for a tally, we use the method ConnectorPlugin::FindMissingVolume() to extract the volume calculated by MC code itself (see note C.1). For MCNP, 2 options are available: the first one is to obtain missing volumes from the MCNP output file1 (the “o” file). Then, if all missing volume have not been found, the second option is used, that is to say a stochastic volume calculation using MCNP itself. For Serpent, only the second option is used (but of course the stochastic volume calculation use Serpent...).

The principle of the stochastic volume calculation for MCNP is to use a spherical shell neutron source geometry directed inwards containing the necessary cells filled with voids. The required volumes are calculated stochastically by MCNP from the number of particle tracks which have intersected them.

In Serpent, no transport is done, just a stochastic sampling (random points over the whole geometry) is performed.

The Serpent way of volume calculation is much faster than the MCNP way, but in general, the number of “source” particles in Serpent should be greater for a same precision (a factor 2 or 4 greater is generally enough). Nevertheless, Serpent is still faster.

The precision on volumes will depend on the “source” neutron number (that can be modified via MURE::SetVolumeNPS()) as well as the relative size of the cell/source neutron sphere. In the usual procedure, a source sphere containing the whole geometry is searched ; However, the search can sometimes fail, in which case the user is required to give the source sphere characteristics manually with the command:

Shape::SetVirtualSphereRadius()

and

Shape::SetVirtualSphereCenter()

It should be noted that it is not necessary for the whole geometry to be included in the “virtual” sphere ; only cells with unknown volumes must be inside this sphere. A volume warning message is given if the relative precision of the worst volume variance is below 0.5% ; this precision may be change with MURE::SetVolumeWarning().

C.1 Important Note

Stochastic volume calculations are used if a tally bin has a volume equal to -1. The tally bin volume is taken from the cell used in this bin. Cell volumes are determined from Shape volumes. Sometimes, we are unable to determine relatively simple Shape volumes because the Included/Disjoint method can’t return a result (e.g., a shape may be included in another one, but Shape::Included() answers 0 (i.e., not included or unknown)). For these cases it is necessary to avoid stochastic volume calculations, whereas if MC code is able to calculate the volume of such cells, we can use the method MURE::GetMCNPVolume(). For a cell without a volume (i.e. -1), if the volume has been calculated and if the cell has no universe number2, then the cell volume is set to what MC code has calculated.

To summarize, the MURE::BuildMCFile() method calls first the MURE::FindMissingVolume() method and afterward the input MCNP file is written.

NOTE: The tally bin volumes should always be provided ; thus, if a group bin or a lattice bin is used, the calculated volume related to that bin...but not to each individual cell composing the bin. Therefore, the individual cell volume remains unknown.

Appendix D
Back-Stage processes: everything you always wanted to know about MURE* (*But were afraid to ask)

D.1 Steady-state power: Tally Normalization

MCNP provides raw tallies normalized by source neutron ; we have, for the sake of simplicity, choose the same normalization method in Serpent.

For reactor applications, in order to simulate a constant reactor power1, one must normalize all tallies by the equivalent number of neutrons per second for a user-given power value. This normalization factor will change as the isotopic composition of the material evolves with time and therefore needs to be calculated at the lowest level of time discretization, i.e., at every sub-step (CRAM/RK ones), and, a fortiori, also every time neutron transport is performed.

  1. Constant power P is defined by the user once and for all throughout the evolution using the gMURE->SetPower() function . This value is translated into the number of fissions per second.

  2. At every MC step, flux ϕ and average reactions cross-sections σiare read from MC output for every nucleus in every cell. Isotopic proportions are read from the initial user-defined composition and for all subsequent steps these values are solved according to their Bateman equations.

  3. At every CRAM/RK sub-step average cross-sections σi, read in the prior MC step, have (so far) been kept constant. New values for the isotopic proportions Ni+1 are solved for using former values of other parameters. P is constant by definition : a new value of the tally normalization factor is therefore calculated in order to keep P constant.

D.1.1 Energy Released from MURE fissions

The energy released per fission is around 200 MeV, but the exact value varies by +/- 5% and varies primarily with the mass of the fissile isotope.

By default MURE assumes that the energy released from fission for each fissionable isotope is different. For several common fissile isotopes, the fission energy values are taken from: W.H. Walerk, “Mass balance estimate of energy per fission in reactors AECL-3109, 1968”. For all other fissionable isotopes, fitted values are used. The linear fit (as a function of mass) is made for the common isotopes [232Th - 241Pu] and extrapolated to the uncommon ones.

For example (see [39]), in the case of 235U, the energy released by fission Q is given in Table D.1.

Table D.1: Energy Released per fission for 235U.
Fission Product kinetic energy 166.2 ± 1.3 MeV
Prompt fission γ 8.0 ± 0.8 MeV
Delayed (after β decays) γ 7.2 ± 1.1 MeV
β decays 7.0 ± 0.3 MeV
Neutron kinetic energy 4.8 ± 0.1 MeV
(n,γ) capture(6 MeV × 1.42 neutrons) 8.5 MeV
Energy released Q 201.7 ± 0.6 MeV

In MURE, the default value is Q235U = 201.7 MeV. The energy of neutrinos (in this case, 9.6 MeV), is not taken into account.

If a user wishes to use specific values, 2 methods are available. The first one, MURE::SetAllFissionEnergies() allows to give a common constant value to each fissionable isotopes (default=200 MeV) ; the second method, MURE::SetFissionReleasedFile() allows user to give a file with the required value for some (or all) fissionable isotopes. Each line of the file is

Z A Qvalue

where Z and A are the proton and nucleon numbers and Qvalue is the energy released by fission (in eV). When using this method, the Q values are taken as explained above (using the linear fit) except for the isotopes specified in this file.

NOTE: MURE::SetAllFissionEnergies() and MURE::SetFissionReleasedFile() must be called BEFORE any Material declaration to be taken into account.

The fact that MURE takes into account different fission energy released values, means that there will be slight differences in the way the normalization factor is calculated, and therefore slight difference in inventory when the fissile nucleus changes as a function of time when compared with the old method.

D.1.2 How the normalization factor is calculated

The total power delivered by the total number of fissions in a given instant (in MC raw units) is given by the double sum shown hereafter :

∑ ∑ PMC = Nijσjiϕj ξ j i MC

where the subscript i refers to the nucleus considered whereas j refers to the cell considered. The sum is to be carried out over all cells and all nuclei belonging to each cell’s material.

PMC can be thought of as the power delivered per source particle transported in the global geometry in question : If ϕj is in standard MC cell tally output units (particles per unit surface per source particle) then PMC will be also normalized per source particle.

Therefore, the scaling factor α, by which tallies will be normalized and which has to be taken into account in order to simulate a constant burn-up at the desired power Puser (if Puseris provided in watts and ξ in eV) is written :

α ≡ ----Puser---- = ∑--∑------Puser---------- PMC 1,6.10−19 i jN jiσjiϕjMC ξ 1,6.10−19

D.1.3 An alternative approach

For a critical system, if one multiplies the desired power in fissions per second times the average number ν of neutrons emitted per fission, one obtains the amount of existing neutrons in the geometry at a given instant. This can be thought of as the intensity of the neutron source needed to achieve this desired power. This value, in neutrons per second, should be equivalent to the value of α calculated in the previous paragraph. If the system is critical (Keff = 1) this number is constant during a time interval small enough so that the amount of fissile material remains unchanged in a first approximation. If the system is critical and Puser is in unit Watt and ξ in unit Electron-Volt, then:

∼ --Puser-- α = ν ξ 1.610−19

For our problem example, the first method gives a value of alpha of roughly 8.28010+12 whereas the second method gives a value of 8.12410+12 with νtaken at 2.599. The two values are therefore comparable within less than 2 percent.

This normalization can be activated by setting

gMURE->GetEvolutionSolver()->SetPowerNormalizationFactorNuOverKeff();

in the input. MURE will add appropriate global tallies to calculate the systems ν.

Bibliography

[1]   O. Méplan, A. Nuttin, O. Laulan, S. David, F. Michel-Sendis, and J. Wilson. MURE: MCNP Utility for Reactor Evolution – Description of the methods, first applications and results. In European Nuclear Society, editor, ENC 2005-European Nuclear Conference. Nuclear Power for the XXIst Century: From basic research to high-tech industry, pages 1–7, 2005.

[2]   F. Michel-Sendis, O. Méplan, S. David, A. Nuttin, A. Bidaud, J. N. Wilson, and O. Laulan. Plutonium incineration and uranium 233 production in thorium fueled light water reactors. In GLOBAL 2005: International Conference on Nuclear Energy Systems for Future Generation and Global Sustainability, page 132, Tsukuba, Japan, October 2005. Tsukuba National Laboratory For High Energy Physics.

[3]   X-5 Monte Carlo Team. MCNP – A General Monte CarloN-Particle Transport Code, Version 5. Technical report, Los Alamos National Laboratory, LA-UR-03-1987, April 2003.

[4]   Jaakko Leppänen, Maria Pusa, Tuomas Viitanen, Ville Valtavirta, and Toni Kaltiaisenaho. The Serpent Monte Carlo code: Status, development and applications in 2013. Annals of Nuclear Energy, 82:142–150, August 2015.

[5]   M. Aufiero, A. Cammi, C. Fiorina, J. Leppänen, L. Luzzi, and M. E. Ricotti. An extended version of the SERPENT-2 code to investigate fuel burn-up and core material evolution of the Molten Salt Fast Reactor. Journal of Nuclear Materials, 441:473–486, 2013.

[6]   J. S. Hendricks, G. W. McKinney, H. R. Trellue, J. W. Durkee, T. L. Roberts, H. W. Egdorf, J. P. Finch, M. L. Fensin, M. R. James, D. B. Pelowitz, and L. S. Waters. MCNPX version 2.6.A . Technical report, Los Alamos National Laboratory report LA-UR-05-8225, 2005.

[7]   D. I. Poston and H. R. Trellue. User‘s manual, version 2.0 for Monteburns, version 5B, Los Alamos National Laboratory report LA-UR-99-4999, 1999.

[8]   L. M. Petrie, W. C. Jordon, and A. L. Edwards. SCALE: a modular code system for performing standardized computer analyses for licensing evaluation. version 5.1 – Volume 1-3, ORNL/TM-2005/39, 2006.

[9]   R. S. Babcock, D. E. Wessol, C. A. Wemple, and S. C. Mason. The MOCUP Interface: A Coupled Monte Carlo/Depletion System. In Topical Meeting on Advances in Reactor Physics, Knoxville, TN, volume III, page 368, April 1994.

[10]   J. Cetnar, W. Gudowski, and J. Wallenius. MCB: A continuous energy Monte Carlo Burnup simulation code. Technical report, Actinide and Fission Product Partitioning and Transmutation, OCDE/NEA, 1999.

[11]   W. Haeck and B. Verboomen. VESTA, An Optimum Approach to Monte Carlo Burnup. Nuclear Science and Engineering, 156(2):180–196, June 2007.

[12]   J. Miss, O. Jacquet, F. Bernard, B. Forestier, W. Haeck, and Y. Richet. First validation of the new continuous energy version of the MORET5 Monte Carlo code. In International Conference on the Physics of Reactors, PHYSOR 2008, 2008.

[13]   Emeric Brun, Eric Dumonteil, and Fausto Malvagi. Systematic uncertainty due to dtatistics in Monte Carlo burnup codes: Application to a simple benchmark with TRIPOLI-4-D. Progress in Nuclear Science and Technology, 2(0):879–885, October 2011.

[14]   Jason Chao. COBRA-3C/RERTR - A Thermal-Hydraulic Subchannel Code with Low Pressure Capabilities. Technical report, Science Applications, Inc., December 1980.

[15]   D. Basile, R. Chierici, M. Beghi, E. Salina, and E. Brega. COBRA-EN, an Updated Version of the COBRA-3C/MIT Code for Thermal-Hydraulic Transient Analysis of Light Water Reactor Fuel Assemblies and Cores. Technical report, OECD/NEA, 1999.

[16]   R. E. MacFarlane and D. W. Muir. The NJOY Nuclear Data Processing System Version 91. Technical report, Los Alamos National Laboratory report LA-12740-M, 1994.

[17]   J. R. Parrington, H. D. Knox, S. L. Breneman, E. M. Baum, and F. Feiner. Chart of the nuclides. fifteenth edition, revised to 1996. Technical report, Knolls Atomic Power Lab., Schenectady, NY (United States), 1996.

[18]   Richard B. Firestone, Coral M. Baglin, and S. Y. Frank Chu, editors. Table of isotopes. Wiley, New York, 1999.

[19]   Mark A. Kellett, Olivier Bersillon, and R. W. Mills. The JEFF-3.1/-3.1.1 radioactive decay data and fission yields sub-libraries. Technical Report NEA 6287, Nuclear Energy Agency, 2009.

[20]   Gaël Guennebaud, Benoît Jacob, et al. Eigen v3. http://eigen.tuxfamily.org, 2010.

[21]   Niels Lohmann. JSON for Modern C++. https://github.com/nlohmann/json, 2013-2021.

[22]   Patrick Boettcher. JSON schema validator for JSON for Modern C++. https://github.com/pboettch/json-schema-validator, 2020.

[23]   American National Standard. Neutron and gamma-ray flux-to-dose rate factors, ANSI/ANS–6.1.1, 1977.

[24]   D. E. Cullen. PREPRO 2019: ENDF/B Pre-Processing Codes. IAEA-NDS-39, Rev. 19.

[25]   A. J. Koning, D. Rochman, S. van der Marck, J. Kopecky, J. Ch. Sublet, S. Pomp, H. Sjostrand, R. Forrest, E. Bauge, H. Henriksson, O. Cabellos, S. Goriely, J. Leppanen, H. Leeb, A. Plompen, and R. Mills. TENDL-2014: TALYS-based evaluated nuclear data library.

[26]   R. E. MacFarlane. The NJOY Nuclear Data Processing System, Version 2016. LA-UR-17-20093.

[27]   M. Oettingen, C. Döderlein, E. D’Agata, K. Tuček, and J. Cetnar. Comparison of MCB and FISPACT burn-up performances using the HELIOS experiment technical specifications. Nuclear Engineering and Design, 242:399–412, January 2012.

[28]   M. J. Bell. ORIGEN: The ORNL Isotope Generation and Depletion Code. ORNL-4628. Oak Ridge National Laboratory, 1973.

[29]   A. G. Croff. ORIGEN2: A versatile computer code for calculating the nuclide compositions and characteristics of nuclear materials. Nucl. Technol., 62:335–352, 1983.

[30]   William Press. Numerical recipes in C : the art of scientific computing. Cambridge University Press, Cambridge Cambridgeshire New York, 1992.

[31]   Maria Pusa and Jaakko LeppÀnen. Computing the Matrix Exponential in Burnup Calculations. Nuclear Science and Engineering, 164(2):140–150, February 2010.

[32]   Maria Pusa. Numerical methods for nuclear fuel burnup calculations. Aalto University, 2013.

[33]   Maria Pusa. Higher-Order Chebyshev Rational Approximation Method and Application to Burnup Equations. Nuclear Science and Engineering, 182(3):297–318, March 2016.

[34]   Maarten Becker. Code and Data Enhancements of the MURE C++ Environment for Monte-Carlo Simulation and Depletion. atw - International Journal for Nuclear Power, 65(6/7):337–340, 2020.

[35]   Timothy A. Davis. Algorithm 832: UMFPACK V4.3—an unsymmetric-pattern multifrontal method. 30(2):196–199.

[36]   E. F. Shores. Data updates for the SOURCES-4A computer code. Nuclear Instruments and Methods in Physics Research Section B: Beam Interactions with Materials and Atoms, 179(1):78–82, June 2001.

[37]   James F. Ziegler and Jochen P. Biersack. The Stopping and Range of Ions in Matter, 1985.

[38]   W. B. Wilson, M. Bozoian, and R. T. Perry. Calculated α-induced thick target neutron yields and spectra, with comparison to measured data. Technical report, 1988.

[39]   M. F. James. Energy released in fission. Journal of Nuclear Energy, 23(9):517–536, November 1969.

1In french, a mûre is a blackberry ; the blue fish with legs of the MURE’s logo is the Darwin evolution symbol (a fish leaving the sea), with a small wrench that symbolize that MURE is a tool. The fish climb on a tree with blackberries ; they represents the nuclei tree need for the evolution, each drupelet of the berries being the nucleons. The small serpent is of course link with the biblical tree of knowledge (but here for the neutronics knowledge) as well as the Serpent Monte-Carlo code used together with MCNP.

2Basic knowledge of C/C++ may help to understand MURE. Nevertheless, careful reading of examples is sufficient to understand and use the package. In appendix 9.3.4.5, a very short introduction to C++ terminology is given in order to facilitate the use of the User Guide and examples.

3For reactor physics, it is not generally necessary to take into account the evolution of every material (such as reflectors, vessels, ...).

4It is, nevertheless, relatively easy to make evolution of user defined MCNP file.

5MURE v1 was designed only for MCNP. It was closely linked to MCNP. The first phase of MURE v2 has been to decouple MURE from MCNP ; the second phase was to couple MURE to a MC code in a “generic” way”. This has caused the break in backward compatibility.

6For csh shell, replace the “export” by a “setenv” and the “=” by a blank char “ “.

7Output files are mainly MCNP/Serpent input files named inp.example_name. All output files have been compressed ; so you have to make a “gunzip *” inside MURE/examples/Static/Reference_Outputs.

8See comments in the file.

9The “r” files of MCNP, in the MCNP version, will be removed (they are in general not useful and very large).

10One can use either, gMURE->SetDATADIR(path) or the environment variable DATADIR (setenv DATADIR path in csh). Note that if you are using gMURE->SetDATADIR(path) you must call it before any Material definition.

11The main difference between B-7.1 and B-6.8 is for the fast fission of Pu-239.

12Most of LINUX distributions provide LAPACK package ; a ’locate liblapack.so’ shows if it is already installed ; if only liblapack.so.version_number is found, you either have to install the “lapack-devel” package or to make a ’ln -s /path_to_lapacklib/liblapack.so.version liblapack.so’ in the MURE/lib directory.

13For some old developers the old Makefile is also available in the MURE/source/src directory.

1In C++ this operator is known as the bitwise right shift operator

2The importance given in the Cell constructor will be the same for all transported particles (define by one of the MURE::SetMode methods) ; if user want to give specific importances for each particle, one needs to call for ALL transported particles the Cell::AddParticle(“particle_Name”, importance_for _that_part).

3i.e. a pointer on a Material, see Appendix 9.3.4.5

4Because we are in a lattice, this whole space is cut by the lattice generator. It is also possible to give a bigger volume to obtain the same kind of result ; however, MCNP has to check all surfaces of this volume for a neutron in the cell so it will result in a longer MCNP run...

1In this figure, we also show the MURE cross-section data base built at IPNO/LPSC ; this base stores nuclei from ENDF, JENDL and JEFF. The relatively small difference in nuclei number between MCNP and MURE bases is due to the fact that this tree has been built at 300K. When choosing higher temperatures (as in standard reactors), the difference will increase because the MURE base has been processed for different temperatures (500K, 600K, 800K, 1000K, 1200K, 1500K) for minor actinides and main fission products. This base is not part of the MURE package.

2This branching ratio is defined as BR = GGSS+m- ; thus the ground state is produced BR ×ReactionRate and isomeric state is produced (1 −BR) ×ReactionRate.

3This branching ratio is defined as BR = -GS-- GS+m ; thus the ground state is produced BR ×ReactionRate and isomeric state is produced (1 −BR) ×ReactionRate.

4One can make ACE binary format from ASCII ACE format using MCNP makxsf exec. Even if NJOY is theoretically able to create such files (ACE binary files), our experience shows that only ASCII files made by NJOY are usable. So you can process ACE ASCII files from ENDF format with NJOY and then, use makxsf to convert them to binary.

5An ASCII file may be used but since the run time increases significantly, we recommend to convert the ASCII ACE file to a binary ACE file with makxsf.

6This threshold has been chosen in order to suppress “exotic” reactions, such as (n,nα), as well as other insignificant reactions and to keep “essential” reactions ; a file, SuppressReaction.dat in the ReactionList directory contains the removed reactions. Don’t forget to verify it!!!!

1On Fedora distributions, the suitesparse-devel package should be installed in order to used this functionality.

2On Fedora or OpenSuSE there is a significant degradation of the performances.

3When using parallel calculation for MCNP (omp,mpi) this factor is reduced: for example, running MCNP in omp on a 8 processor machine leads to only 16 times faster for a simple problem.

4This is of the same order of magnitude than the deviation observed in a series of N identical evolution starting with a different random seed. Using this multi-group method, unresolved probability tables are not taken into account for reaction rates evaluations ; this explain that for fast systems, the discrepancy is larger than for thermal ones.

5The default is 4, and it can be changed with MURE::SetFitRangeNumber().

1But if FP and/or AM check boxes are selected, they are not part of the sum.

2The input cooling time is then chosen as the lower bound of the time binning defined in the “Evolution period” frame ; see the exact value in the terminal window.

1Except if the ReactorAssembly is used to fill a core lattice.

1One cannot write void AnOtherMethod(double y=2, double x).

1The principle of this method is to “run” MCNP in “information mode”, then search the print table 60 containing the cell volumes information.

2The volume given in table 60 is the volume of the cell alone ; if a universe is defined, then the volume of the true cell is the volume of all places where that universe is input. Thus only stochastic volume calculations can be made.

1For this purpose, only fissions are considered as source of power, at an average value of 200 MeV per fission.