NOrmalizing Tally to get Flux value

User Support
17

Well, I guess that there will be still someone who has these problems about the normalization of neutron flux.

Actually, Paul and other friendly guys have given many detailed explanations, and I could finally solve the problem because of these useful suggestions :grin:

And I noticed that there are many posts about the normalization of neutron flux, so, I’ll give an example to show how to normalize the neutron flux under depletion for a k-eigenvalue calculation, this example is based on Paul’s example Tally a flux spectrum and OpenMC’s example Pincell Depletion.

Ok, let’s begin :wink:

Define materials

#Import some necessary packages.
%matplotlib inline
import openmc
import math
import numpy as np
from matplotlib import pyplot

# Define the fuel
uo2 = openmc.Material(name="uo2")
uo2.add_element("U", 1, percent_type="ao", enrichment=4.25)
uo2.add_element("O", 2, percent_type="ao")
uo2.set_density("g/cm3", 10.4)
uo2.temperature = 670 # 400 Celsisus

# Define the helium gap
he = openmc.Material(name='gap')
he.add_element('He', 1.0)
he.set_density('g/cm3', 0.0014)

# Define the zirconium cladding
zr4 = openmc.Material(name="clad")
zr4.add_element('Zr', 0.982, percent_type="wo")
zr4.add_element('Sn', 0.015, percent_type="wo")
zr4.add_element('Fe', 0.002, percent_type="wo")
zr4.add_element('Cr', 0.001, percent_type="wo")
zr4.set_density("g/cm3", 6.5)
zr4.temperature = 533 # 260 Celsisus

# Define the coolant
water = openmc.Material(name="water")
water.add_element("O", 1)
water.add_element("H", 2)
water.set_density("g/cm3", 1.0)
water.add_s_alpha_beta("c_H_in_H2O")
water.temperature = 533 # 260 Celsisus 

materials = openmc.Materials([uo2, he, zr4, water])

Define the geometry

# Basic geometry parameters
pitch = 1.3 # unit is cm
fuel_outer_radius = openmc.ZCylinder(r=0.4) 
clad_inner_radius = openmc.ZCylinder(r=0.43)
clad_outer_radius = openmc.ZCylinder(r=0.5)
left = openmc.XPlane(x0=-pitch/2, boundary_type='reflective')
right = openmc.XPlane(x0=pitch/2, boundary_type='reflective')
bottom = openmc.YPlane(y0=-pitch/2, boundary_type='reflective')
top = openmc.YPlane(y0=pitch/2, boundary_type='reflective')

# Regions
fuel_region = -fuel_outer_radius
gap_region = +fuel_outer_radius & -clad_inner_radius
clad_region = +clad_inner_radius & -clad_outer_radius
water_region = +left & -right & +bottom & -top & +clad_outer_radius

# Fill the regions with materials
fuel = openmc.Cell(name='fuel')
fuel.fill = uo2
fuel.region = fuel_region
gap = openmc.Cell(name='helium gap')
gap.fill = he
gap.region = gap_region
clad = openmc.Cell(name='clad')
clad.fill = zr4
clad.region = clad_region
moderator = openmc.Cell(name='moderator')
moderator.fill = water
moderator.region = water_region

# Universe
root_universe = openmc.Universe(cells=(fuel, gap, clad, moderator))

geometry = openmc.Geometry()
geometry.root_universe = root_universe

root_universe.plot(width = (pitch, pitch))

Define the setting

settings = openmc.Settings()

settings.temperature = {'method': 'interpolation'} 
settings.particles = 2000
settings.generations_per_batch = 10
settings.inactive = 50
settings.batches = 200

uo2.volume = math.pi * 0.4 ** 2

# Finally, don't forget to export to xml.
materials.export_to_xml()
geometry.export_to_xml()
settings.export_to_xml()

Tally neutron flux with energy filter

# Create equal-lethargy energies to put in filter
energies = np.logspace(np.log10(1e-5), np.log10(20.0e6), 501)
e_filter = openmc.EnergyFilter(energies)

# Create tally with energy filter
tallies = openmc.Tallies()

tally = openmc.Tally(name='flux')
tally.filters = [e_filter]
tally.scores = ['flux']
tallies.append(tally)

# No filters here
fiss_rate = openmc.Tally(name='fiss_rate')
neu_fission = openmc.Tally(name='neu_fission')
fiss_rate.scores = ['fission']                   # Total fission reaction rate [fission/source]
neu_fission.scores = ['nu-fission']              # Total production of neutrons due to fission [neutrons/source]
tallies += (fiss_rate, neu_fission)

# Export to "tallies.xml"
tallies.export_to_xml()

Setting up for the depletion

import openmc.deplete

chain = openmc.deplete.Chain.from_xml("./chain_casl_pwr.xml")  # This library is simplified.
operator = openmc.deplete.Operator(geometry, settings, "./chain_casl_pwr.xml")

power = 174  # [W/cm]
time_steps = [30 * 24 * 60 * 60] * 6  # These stpes are not cumulative.

integrator = openmc.deplete.PredictorIntegrator(operator, time_steps, power)
integrator.integrate()

Post Processing

After the running, we will get 7 files, one for each our of 6 depletion steps and the initial state. Pranto said that “Depletion module will write tally data to statepoint files for each burnup steps in openmc_simulation_n.h5, current time step.By opening each statepoint file, you will get neutron flux for different burnup step.”

So, let’s normilize the neutron flux in the end of life step.

# To load data from the statepoint file in the EOL burnup step
sp = openmc.StatePoint('openmc_simulation_n6.h5')

# Get the mean values for each energy bin by using the mean attribute on the tally.
t = sp.get_tally(scores=['flux'])
flux = t.mean.ravel()        # [neutrons-cm/source]
flux_unc = t.std_dev.ravel()

# The unit of flux is [n/cm-src], we need to convert it into [n/cm2-s]
fiss = sp.get_tally(name='fiss_rate')
fiss_rate = fiss.mean.ravel()
neu = sp.get_tally(name='neu_fission')
neu_fission = neu.mean.ravel()
V = uo2.volume               # [cm^3]
P = 59.5                     # [J/sec], actually, I'm not sure.
Q = 3.2e-11                  # [J/fission]
nu = neu_fission/fiss_rate   # [neutrons/fission]
k = neu_fission              # [neutrons/source]
Flux = flux/V*(P*nu/(Q*k))  # [neutrons/cm^2-sec]

# We can convert some variables into matlab variables.
sio.savemat('Flux_EOL.mat', {'Flux': Flux})
sio.savemat('Energy.mat', {'energies': energies})

# Obtain the width of Lethargy Energy
num_groups=len(Flux)
delta_LethEnergy=np.zeros(num_groups)
normalized_flux=np.zeros(num_groups)

for i in range(num_groups):
    delta_LethEnergy[i]=math.log(energies[i+1]/energies[i])

# Calc normlized flux
integral_flux=0

for i in range(num_groups):
    integral_flux=integral_flux+Flux[i]

for i in range(num_groups):
    normalized_flux[i] = Flux[i]/integral_flux/delta_LethEnergy[i]

#  Use matplotlib to plot the normalized flux versus the energy.
pyplot.semilogx(energies[:-1],normalized_flux)
pyplot.grid()
pyplot.xlabel('Energy [eV]')
pyplot.ylabel('Normalized flux/lethargy width')

Matters needing attention

The unit of flux in the tally is [neutron-cm/source], we need to convert it to [neutron/cm^2-sec]

  • flux — [neutron-cm/source]

  • flux/V = [neutron/cm^2-source]

    P – reactor power of your element [J/s]

Q = 200 MeV for U-235 or 3.2x10-11 [J/fission]

nu – score ‘nu-fission’/'fission reaction rate' [neutrons/fission]

k – eigenvalue [neutrons/source]

V — volume of your element/target

  • normalization factor = Pnu/ (Qk) = [J/sec * neutrons/fission / (J/fission * neutrons/source) ] = [source/sec]

  • Normalized Flux = flux * normalization factor = [neutron/cm^2-sec]

  • Also, I gave the code about how to get the Normalized flux/lethargy width.

Some detailed discussions or documents about this topic.

  1. OpenMC’s Document about Normalization of tally results
  2. Question about power designation while depleting
  3. OpenMC’s Normalization Factor (This one gives an detailed example)
  4. Normalization Factor in Section 8.3
  5. Unit transformation in a fixed source problem
  6. To find neutron flux, atomic density, and burn-up level through depletion calculation

Finally

There may be some mistakes about the detailed explanation that I’m not sure, but don’t hesitate to put them forward. :yum:

15 Likes