<figure><img src="image_40.png" /><figcaption>Figure 8.10: MureGui Radiotoxicity Tab</figcaption></figure>

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

This file is saved in the *tox/DATA_000*. The activity of this initial composition is obtained by doing

MureGui -type a tox

One can see that the activity starts from 1 mCi and reaches secular equilibrium (14 mCi) after ~ 2 × 10⁶ years.

**Remarks:**

* The "type -a" is necessary because the *DATA_000* file is in *ASCII* and not in binary (the default of *MureGui*)

* Depending on your installation, you probably need to set the environment variable *DATADIR* to the correct path (for example in *csh*, *setenv DATADIR /MURE_install_path/MURE/data*)

### 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).

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 \pm \Delta E$ enters in combo boxes ($E$ and $\Delta E$ are both in keV, see fig. 8.11).

**“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)

* ASCII format: dump the last plotted data in a file.

* “Source Input for MURE” format: a *MCNPSource C++* form file is generated in order to be used for a *MURE/MCNP* simulation. This source as the energy spectra of the selected plot, it is a collimated, uniformly distributed on a disk of a radius of 10cm source. Particles are emitted perpendicularly to the disk. This can be very useful for radio-protection studies.

**“Use binning from file” button:** This button is used to load the energy binning for Spectrum from a file. The format of this file is the same of the one used for Spectrum class (see for example `/path_to_MURE/MURE/gui/NEA_Binning_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 [m³]
```

Finally open it with *MureGui* :

Write in your 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*, ...)
<figure><img src="image_41.png" /><figcaption>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, ...).</figcaption></figure>
### 8.4.1 Spectrum Class

It's a base class to define energy spectra for $\gamma$, $\beta$, $\alpha$ and neutrons. This class defines Spectrum objects, i.e. one array store intensities of other store energies. The size of these array are `NumberOfBins+1`. To declare a Spectrum proceed as follow:

```cpp
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:

```cpp
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:

```cpp
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 another with `Spectrum::AddSpectrum()` but **only if the two spectra have the same energy binning**:

```cpp
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 $\gamma$).

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):

```cpp
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" ($\alpha$ decay), "B" ($\beta^-$ decay), "E" (electronic capture) and "I" (Isomeric Transition)).

`GammaSpectrum::ReadENSDF()` searches in `ENSDF` files the gamma transition of the $A_ZX$ daughter(s) from these values: decay constant of $A_ZX$ ($\lambda$), 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 = \lambda \times N_{A_Z X} \times RI \times NR \times BR/100 \ [bq]$$
where $RI \times NR$ is the number of photon per 100 decays of the parent for this decay branch. $RI \times NR \times 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 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)

$$ \text{EndpointEnergy} = \text{Qvalue} + \text{EnergyOfTheDecayingLevel} - \text{EnergyOfTheLevelOfTheDaughter} $$

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) = \text{Atoms} \times \lambda \times RI \times NR \times BR \times f(E) $$

where

$$ f(E) = \text{Norm} \times \frac{2\Pi\eta}{1 - \exp(-2\Pi\eta)} \times W \times \sqrt{W^2 - 1} \times (\text{EndpointEnergy} - E)^2 $$

$$ \text{with } \eta = \frac{Z_{\text{daughter}} \times e^2}{4\Pi\varepsilon_0 h V_e}, \quad V_e = c \times \sqrt{1 - \frac{1}{(1+W)^2}} \quad \text{and} \quad W = 1 + \frac{E}{m_e c^2}. $$

### 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 $\beta^{-}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) = \text{Norm. exp}(-\frac{E}{a}) \cdot \sinh(\sqrt{b \cdot E}) \text{ where } a \text{ and } b \text{ are coefficient depending on the nucleus.} $$

Finally the neutron activity of the nucleus $A_Z X$ by spontaneous fission at energy $E$ is:

$$ A(E) = \lambda \times N_{A_X} \times \langle \nu \rangle \times BR \times W(E) $$

where $BR$ is the branching ratio of the spontaneous fission and $\langle \nu \rangle$ the average neutron per fission. The parameters $\langle \nu \rangle$, $BR$, $a$ and $b$ are available in literature and those used in MURE are located in `/PathToMure/MURE/data/SFNeutronSpectrum/` [24].
### 8.4.5.2 Neutron from $(\alpha, n)$ reactions

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

The stopping power was calculated for a $\text{UO}_2$ media using SRIM[23], $(\alpha, n)$ cross section are from JENDL (JENDL (α,n) Reaction Data File 2005 (JENDL/AN-2005)). The probability $p_i$ of an alpha of energy $E_\alpha$ having a $(\alpha, n)$ reactions on target $i$ while alpha slowing down on $dx$ is $p_i = N_i\sigma_i dx$, where $\sigma_i$ is the $(\alpha, n)$ cross section on target $i$. The probability $p_i$ can be expressed as a function of the variation of energy of the alpha particle on $dx : p_i = N\sigma\frac{dx}{dE}dE$. By integrating the energy of the alpha particle from its initial energy to 0:

$$P_i(E_\alpha) = N_i \int_0^{E_\alpha} \frac{\sigma_i(E)}{-\frac{dE}{dx}} dE$$

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 [21] for more information).

### 8.4.5.3 Neutron from $(\beta^{-}n)$ decays

* This delayed neutron production has been implemented only for $(\beta^{-}n)$ decays; $(\beta^{-}2n)$ decays are not yet taken into account.

* The neutron production is taken into account only when both $(\beta^{-}n)$ decays and neutron energy $E_n$ is given in ENSDF files. In July 2012, only 6 identified nuclei have this two information: $^{85}\text{As}$, $^{87}\text{Br}$, $^{88}\text{Br}$, $^{95}\text{Rb}$, $^{135}\text{Sb}$ and $^{137}\text{I}$.

## 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).

* First define the sample Material than can emits source particles $(\gamma, n, \beta)$.

* Then define the Spectrum (GammaSpectrum, NeutronSpectrum or BetaSpectrum) with one of the constructors

* Finally proceed as follow (for a punctual and isotropic source):

```
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.

* Geometry builders: *GenericReactorAssembly* (abstract class), and the 2 concrete implementations *MCNP::ReactorAssembly* and *Serpent::ReactorAssembly*

* Thermal-hydraulics flow simulation class: *COBRA_EN* class.

* **BATH** (Basic Approach of Thermal Hydraulics): *ThermalCoupling*, *Thermal reader*.

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 cuboid 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**[12, 13] (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:
* only one guide tube pin type is allowed (same dimension for all guide tubes); a guide tube pin must be composed by an inner coolant and a clad layer and a surrounding (coolant) material.

* only one control rod pin type; a control rod pin must be composed by 4 layers (the rod, its clad, a coolant ring ("gap") and outer clad (the guide tube clad)) and a surrounding (coolant) material. The control rod outer clad should be the same as the guide tube pin clad.

* 2 types of fuel pins are allowed: standard fuel pins (a fuel inner part, its clad and the surrounding (coolant) material) or special fuel pins with 3 layers (fuel inner part, a ring of non fissile material (Gd, ...), and its clad) with a surrounding (coolant) material.

**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:

* heat conduction is calculated in homogeneous isotropic media (temperature only depends on radius)

- simulations are done on average elementary channel(s) (no cross-flow). Several "average" channels can be simulated but thermal-hydraulics interactions between them cannot be calculated.

- The code only simulates a liquid phase (no boiling).

## 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 as 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 a transient analysis.

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

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).

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.
<figure><img src="image_42.png" /><figcaption>Figure 9.1: Meshing grids.</figcaption></figure>

<figure><img src="image_43.png" /><figcaption>Figure 9.2: 4 Radial zones and 10 axial levels.</figcaption></figure>
### 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 provided in `RA_1ass.cxx` (1 hexagonal assembly), `RA_1ass_septent.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*, int, int)`

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 called<fnref n="1" /> 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
```

<footnote n="1">¹Except if the `ReactorAssembly` is used to fill a core lattice.</footnote>
If equal size upper and lower plenums are needed, user should call `ReactorAssembly::SetPlenum(thickness, Plenum-Material)` 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/\sqrt{2}$ for a cuboid assembly of side $S$ and to $2 \times 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.

### 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:
<figure><img src="image_44.png" /><figcaption>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.</figcaption></figure>
<figure><img src="image_45.png" /><figcaption>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).</figcaption></figure>

```cpp
...
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).

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

### 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)
PWR 17x17 Assembly

<figure><img src="image_46.png" /><figcaption>Figure 9.5: Example of an assembly created with ReactorAssembly.</figcaption></figure>
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 around 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 these 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.

* `INFILE` is the COBRA-EN input file

* `OUTFILE` is the COBRA-EN output file

* `Mesh_MURE_COBRA_Data.txt`: are written the lattice positions for each pin in the mesh grid of neutronic code and thermalhydraulics code

* `MCNPCellsTemp.data`: are written averages temperatures in each zone and level for fuels, cladding and coolant (K)

* `PowerDeposits.data`: are written the local power deposits and its absolute error for a representative fuel element in each zone and level (total power released in the zone/level is normalized by the number of fuel pin in the zone - unit: W)

* `PLOTFILE` is the COBRA-EN output plotfile

* A file called `yyy_MeshZoneData.txt` is also created in the initial running directory: the number of fuel pins, guide tubes, special pins for each zone are quoted; an equivalent of the mesh `MCNP_Zones` is also written.
<figure><img src="image_47.png" /><figcaption>Figure 9.6: Geometry visualization</figcaption></figure>

## 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/m²) :

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/C4Assembly (c.f. figure 9.6) and some results :

* Power distribution: figure 9.7a

* Exit coolant temperature profile: figure 9.7b.
<figure><img src="image_48.png" /><figcaption>Figure 9.7: Complex COBRA coupled simulation : 20 axial levels and 4*7 radial regions</figcaption></figure>

### 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:

* the capability to be coupled with geometry that has already been built

* it is an explicit way of coupling, despite the fact it takes longer, users can completely control what is being simulated (it is not a "black box")

* capability to simulate a lot of coolant, cladding or fuels types: users just need to add them (thermal data and equations).

#### 9.3.4.1 Capabilities

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

* as fuel: uranium oxide and plutonium/uranium mixed oxide and thorium oxide (full or annular pellet)

* as cladding: zircaloy 4 or 2 and stainless steel 316 SS

* as coolant: light water, heavy water and sodium.

#### 9.3.4.2 What is solved

* heat and mass transfer:

$$ \dot{m}C_p\Delta T = q_pA $$

* steady-state forced convection in case of turbulent flows:

$$ T_p - \bar{T}_{\text{coolant}} = \frac{q_p}{h} $$

* steady-state conduction inside the cladding and the fuel: Fourier law
* radiation exchange between cladding and fuel surfaces is implemented but not used (due to a leak of precision on gas composition and evolution of the gas space dimension during burn-up): *Stefan law*

* indication of pressure losses (due to acceleration, gravity and friction - loss due to peculiarity like grids are not calculated because it is considered as empirical data): acceleration is zero because the velocity is assumed constant, only gravity and friction (Darcy frictions model) are calculated

* thermal conductivity of zircaloy and steel, data on water, sodium, UOX and MOX are from literature.

* Nusselt dimensionless number calculation:

- water: Dittus-Boelter equation $Nu = 0.023 Re^{0.8} Pr^{0.4}$

- sodium: Notter & Sleicher $Nu = 6.3 + 0.0167 Re^{0.85} Pr^{0.93}$

### 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 :

* on materials: flag *SetTemperatureEvolution()* and only on fuel material *SetEvolution()*

* on cells: flag(s) *SetTHLevelPosition()* (the first one at the bottom and then increasing height) and *SetTH-ZonePosition()* if more than one average channel is needed.

**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).
<figure><img src="image_49.png" /><figcaption>Figure 9.8: BATH output example</figcaption></figure>
# 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

```cpp
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:

```cpp
#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`:

```cpp
class Point
{
```
```cpp
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

```cpp
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:

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

P can be translated of $\vec{V}(1, 1)$:

```cpp
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:

```cpp
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

```cpp
Point P(0,4);
```

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

```cpp
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

```c
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 argument<fnref n="1" />. If one writes

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

then a call to

```c
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:

```c
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.

<footnote n="1">¹One cannot write `void AnOtherMethod(double y=2, double x)`.</footnote>
## A.5 Inheritance

One very powerful tools of classes is the *inheritance* mechanism. 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:

```cpp
#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;}
        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)`).

```cpp
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

```cpp
#include <iosstream>
int main()
{
    ...
    std::cout << 'bla bla bla' << std::endl;
    ...
}
```

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

```cpp
#include <iosstream>
using namespace std;
int main()
{
    ...
}
```
```cpp
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",

```cpp
#include <iosstream>
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.

<figure><img src="image_50.png" /><figcaption>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).</figcaption></figure>
<figure><img src="image_51.png" /><figcaption>Figure B.2: The tree of figure B.1 after simplification. Since one node was included in the other, it has been removed.</figcaption></figure>

<figure><img src="image_52.png" /><figcaption>Figure B.3: An other example of tree before simplification.</figcaption></figure>
<figure><img src="image_53.png" /><figcaption>Figure B.4: The tree of figure B.3 after simplification.</figcaption></figure>
# 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 *ConnectorPlu-gin::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 file<fnref n="1" /> (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 changed with *MURE::SetVolumeWarning()*.

<footnote n="1">¹The principle of this method is to “run” MCNP in “information mode”, then search the print table 60 containing the cell volumes information.</footnote>
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 number<fnref n="2" />, 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.

<footnote n="2">²The 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.</footnote>
# 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 power<fnref n="1" />, 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 RK step, 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 $\phi$ and average reactions cross-sections $\sigma_i$ are 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 RK step average cross-sections $\sigma_i$, read in the prior MC step, have (so far) been kept constant. New values for the isotopic proportions $N_{i+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.

<footnote n="1">¹For this purpose, only fissions are considered as source of power, at an average value of 200 MeV per fission.</footnote>
<table><caption>Table D.1: Energy Released per fission for $^{235}$U.</caption><tr><td>Fission Product kinetic energy</td><td>166.2 ± 1.3 MeV</td></tr><tr><td>Prompt fission γ</td><td>8.0 ± 0.8 MeV</td></tr><tr><td>Delayed (after β decays) γ</td><td>7.2 ± 1.1 MeV</td></tr><tr><td>β decays</td><td>7.0 ± 0.3 MeV</td></tr><tr><td>Neutron kinetic energy</td><td>4.8 ± 0.1 MeV</td></tr><tr><td>(n, γ) capture(6 MeV × 1.42 neutrons)</td><td>8.5 MeV</td></tr><tr><td>Energy released Q</td><td>201.7 ± 0.6 MeV</td></tr></table>

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 [20]), in the case of $^{235}$U, the energy released by fission Q is given in Table D.1. In MURE, the default value is $Q_{235U} = 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 :

$$P_{MC} = \sum_{j} \sum_{i} N_{i}^{j} \sigma_{i}^{j} \phi_{MC}^{j} \xi$$

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.

$P_{MC}$ can be thought of as the power delivered per source particle transported in the global geometry in question : If $\phi^j$ is in standard MC cell tally output units (*particles per unit surface per source particle*) then $P_{MC}$ will be also normalized per source particle.

* $N_i^j$ is the number of nuclei *i* present in cell *j*.
* $\sigma_i^j$ is the average fission cross-section of nucleus i in cell j.

* $\phi^j$ is MC's neutron flux track-length estimation for cell j.

* $\xi$ is the average energy delivered in one fission (~200 MeV).

Therefore, the scaling factor $\alpha$, 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 $P_{user}$ (if $P_{user}$ is provided in watts and $\xi$ in eV) is written :

$$ \alpha \equiv \frac{P_{user}}{P_{MC} \ 1,6.10^{-19}} = \frac{P_{user}}{\sum_{i} \sum_{j} N_{i}^{j} \sigma_{i}^{j} \phi_{MC}^{j} \xi \ 1,6.10^{-19}} $$

### D.1.3 An alternative approach (for information only)

For a critical system, if one multiplies the desired power in fissions per second times the average number $\nu$ 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 $\alpha$ calculated in the previous paragraph. If the system is critical ($K_{eff} = 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 $P_{user}$ is in watts and $\xi$ in electron-volts, then:

$$ \alpha \approx \nu \frac{P_{user}}{\xi \ 1.6 \ 10^{-19}} $$

For our problem example, the first method gives a value of alpha of roughly $8.280 \ 10^{+12}$ whereas the second method gives a value of $8.124 \ 10^{+12}$ with $\nu$ taken at 2.599. The two values are therefore comparable within less than 2 percent.
Bibliography

[1] Méplan O., Nuttin A., Laulan O., David S., Michel-Sendis F. et al., “MURE : MCNP Utility for Reactor Evolution - Description of the methods, first applications and results”, Proceedings of the ENC 2005 (CD-Rom) - ENC 2005 - European Nuclear Conference. Nuclear Power for the XXIst Century : From basic research to high-tech industry, France

[2] Michel-Sendis F., Méplan O., David S., Nuttin A., Bidaud A. et al., “Plutonium incineration and uranium 233 production in thorium fueled light water reactors”, GLOBAL 2005 Proceedings (CD-Rom) - GLOBAL 2005: International Conference on Nuclear Energy Systems for Future Generation and Global Sustainability, Japan

[3] X-5 Monte Carlo Team, “MCNP — A General Monte Carlo N-Particle Transport Code, Version 5”, Los Alamos National Laboratory report LA-UR-03-1987 (April 2003)

[4] SERPENT, PSG2 / Serpent Monte Carlo Reactor Physics Burnup Calculation Code, 2011.
<http://montecarlo.vtt.fi>

[5] Aufiero, M., Cammi, A., Fiorina, C., Leppänen, J., Luzzi, L., and Ricotti, M. (2013a) "An extended version of the Serpent-2 code to investigate fuel burn-up and core material evolution of the molten salt fast reactor." J. Nucl. Mat., 441 (2013) 473-486.

[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”, Los Alamos National Laboratory report LA-UR-05-8225 (2005)

[7] H.R. Trellue and D.I. Poston, “User’s Manual, version 2.0 for MONTEBURNS, version 5B”, Los Alamos National Laboratory report LA-UR-99-4999 (1999)

[8] “SCALE: A Modular Code System for Performing Standardized Computer Analyses for Licensing Evaluations”, ORNL/TM-2005/39, Version 5.1, Vols. I-III, (2006).

[9] R. S. Babcock, D. E. Wessol, C. A. Wemple, and S. C. Mason, “The MOCUP Interface: A Coupled Monte Carlo/Depletion System”, 1994 Topical Meeting on Advances in Reactor Physics, Knoxville, TN, p. III-368 (April 11-14, 1994)

[10] J. Cetnar, W. Gudowski and J. Wallenius, “MCB: A continuous energy Monte Carlo Burnup simulation code”, In “Actinide and Fission Product Partitioning and Transmutation”, EUR 18898 EN, OECD/NEA (1999) 523.

[11] R. E. MacFarlane and D. W. Muir, “The NJOY Nuclear Data Processing System Version 91”, Los Alamos National Laboratory report LA-12740-M, (October 1994).
[12] Jason Chao, "COBRA-3C/RERTR - A Thermal-Hydraulic Subchannel Code with Low Pressure Capabilities", Science Applications, Inc. (December 25, 1980)

[13] D. Basile, R. Chierici, M. Beghi, E. Salina and E. BregaCOBRA-EN, an Updated Version of the COBRA-3C/MIT Code for Thermal-Hydraulic Transient Analysis of Light Water Reactor Fuel Assemblies and Cores Report 1010/1 (revised 1.9.99)

[14] J.R. Parrington et al., "Nuclides and Isotopes", 15th edition, General Electric Nuclear Energy, 1996

[15] R.B.firestone, "Table of Isotopes", 8th edition, V.S.Shirley editor, John Wiley & Sons, Inc, 1996

[16] M.A. Kellett, O. Bersillon, R.W. Mills, "The JEFF-3.1/-3.1.1 radioactive decay data and fission yields sub-libraries", NUCLEAR ENERGY AGENCY, NEA No. 6287, OECD 2009, ISBN 978-92-64-99087-6

[17] W. Haeck, B. Verboomen, "An optimum Approach to Monte Carlo Burn-Up", Nucl. Sci. Eng., 156, pp 180-196 (2007)

[18] J. Miss, O. Jacquet, F. Bernard, B. Forestier, W. Haeck, Y. Richet, "First validation of the new continuous energy version of the MORET5 Monte Carlo code", PHYSOR 2008

[19] E. Brun, E. Dumonteil and F. Malvagi, "Systematic Uncertainty Due to Statistics in Monte Carlo Burnup Codes: Application to a Simple Benchmark with TRIPOLI - 4 - D", Progress in Nuclear Science and Technology, Vol. 2, pp.879-885 (2011)

[20] M.F. James, Journal of Nuclear Energy, 23, 517 (1969)

[21] W.B. Wilson, M. Bozoian and R.T. Perry (1988) "Calculated α induced thick target neutron yields and spectra, with comparison to measured data"

[22] Neutron and Gamma-Ray Fluence to Dose Factors, American National Society, ANSI/ANS-6.1.1-1977

[23] James F. Ziegler, "The Stopping and Range of Ions in Matter (SRIM)", http:www.srim.org

[24] E.F. Shores, "Data update for SOURCE-4A computer code", Nuc. Inst. Meth. B, 179, (2001), pp 78-82

[25] M. Oettingen et al., "Comparison of MCB and FISPACT burn-up performances using the HELIOS experiment technical specifications", Nuc. Eng. and Design, 242, (2012), pp 399-412