NOrmalizing Tally to get Flux value

User Support
1

Hye everybody,

i have performed a calculation for IAEA benchmark reactor and i need to get flux in various direction like axial profile and radial profiles. The tally i used in tally.xml is shown below.

1 1 30 -0.5 -0.5 -35 0.5 0.5 35 flux

for example if i get a value of 0.225861832063,0.036718304271,0.053172569348 for thermal, epitehrmal and fast flux for above mesh setting. Can anyone guide me how to normalize these values to get flux values? Thanks to everyone and sorry if this answer to this question is too obvious.

2

The tally results are integrals over the entire phase space - that includes the volume of the mesh. So if you want a value of flux in [n/(cm^2-s)] then you will need to divide by the volume of each mesh cell. Is that what you were looking for?
Adam

3

More or less. In MCNP we multiply F4 tally with power and nu and divide them by 1.6022E-13 x 193 x keff. This gives neutron per cm2 per sec. But here how do we get units of flux balanced

4

Maybe I’m not understanding correctly what you want, but to my knowledge the OpenMC flux score outputs the same units as the MCNP F4 tally. So if what you do for MCNP suffices, then it should here as well.

Is my understanding correct: you want to take the OpenMC flux score and take it from the normalized values that OpenMC outputs in to the mesh element averaged values of the flux that should be expected when operating at power P, right?

5

thank you nelson for answering. Maybe i am confused by your answer in another post on the forum. there you wrote and i quote

“Hey Anthony,
The tally responses are integrated over volume, and so are [# cm^3]…”

so if the units of tally are #-cm^3 then multiplying by power and nu and dividing them by 1.6022E-13x 193 x keff x volume of mesh, it would give #/sec i guess. so my question is do we have to divide the calculation done above by area or not to get the required units of (#/cm^2-sec)?

6

The tally is not normalized by volume, so the units are something like [neutrons-cm/source]. This can be understood quite easily: in a collision estimator tally, OpenMC accumulates 1/Sigma_t which has units of cm, and similarly for a tracklength tally, OpenMC accumulates the tracklength in cm. When you divide by volume, you get [neutrons/cm^2-source]. The normalization factor Pnu/(Qk) gives you [J/secneutrons/fission/(J/fissionneutrons/source)] = [source/sec] so that when you multiply the flux by the normalization factor, you get [neutrons/cm^2-sec].

1 Like
7

thank you very much Paul :slight_smile:

8

Thank you,Paul. But how can I get P with a unit of (J/sec) ? I also notice that we can get a unit of fission/sec by P/Q, is this what we get by fission tally ?

在 2014年8月27日星期三 UTC+8上午1:49:12，Paul Romano写道：

9

P is the reactor power and is generally known a priori. For example, if you have a 1000 MWe reactor with a thermal efficiency of 33%, the thermal power is 3000 MW = 3e9 J/s. P/Q would tell you how many fission reactions per second are necessary to sustain the given power. Once you account for nu (neutrons/fission), then you have the number of neutrons needed to sustain a given power level. All tallies in OpenMC are ‘per source neutron’, so you need to determine how many neutrons/sec there are in your system of interest using a combination of these factors.

10

Hi Paul,

I am interested in expressing the flux in “neutrons/cm2s”. I was reading this post where it is explained how to do it. I would like to know if I understood the normalization factor that you presented [Pnu/(Qk)]:

P – reactor power [J/s]

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

nu – score ‘nu-fission’ [neutrons/fission]

k – eigenvalue [neutrons/source]

Is my interpretation okay?

Thanks,

Javier

1 Like
11

All that is correct except for nu. While the units are indeed [neutrons/fission], the ‘nu-fission’ score doesn’t quite give you what you need. ‘nu-fission’ gives you the total fission neutron production, so you need to divide by the fission reaction rate in order to get the average number of neutrons per fission.

Best,
Paul

4 Likes
12

Thanks Paul.

Javier

13

Paul,

Continuing with this post, a new question comes to my mind.

If I am modeling one fuel element, in the normalization factor, is “P” the total reactor power or the total reactor power divided by the number of fuel elements?

Thanks,
Javier

14

If modeling a single fuel element, then you would use the power of just one element as opposed to the power for the whole core.

Best,
Paul

15

Hello, Paul,

I notice that in a flux spectrum, the ordinate is usually Normalized flux/lethargy width .

So, after we get the flux in [neutrons/cm^2-sec], how do we get the Normalized flux/lethargy width ?

Thank you,
Deng.

16

Sorry to bother, Paul.

I’ve solved it.

Thanks,

Deng.

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:

14 Likes
18

Hi everyone,

As a follow-up on this thread. If I am comparing the flux tally result between a criticality run and a fixed-source run for the same core configuration (assumed subcritical) and I want to retain the final result to be per source particle. Are they directly comparable? Or the tally result of the eigenvalue problem should be multiplied by keff (fission neutron / source neutron) since - I assume similar to MCNP - the result is actually per fission neutron?

Thank you in advance.

19

I’m not sure you could really get tally results that are exactly equivalent between the two. If you are running a fixed-source problem with subcritical multiplication, the flux will be amplified by a factor of ~1/(1-keff). This is contrast to a k-eigenvalue calculation where production of fission neutrons is deliberately scaled by 1/keff to obtain a balanced system.

20

Hi all,

What is the difference between normalizing the tally with Pnu/(Qk) and heating score (from Section 8.3)?

I calculated fission rate with both methods but somehow the latter one did not give me the same total power from what I have set.

Would appreciate any comments!
Kristina

1 Like