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
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 keigenvalue calculation, this example is based on Paulâs example Tally a flux spectrum and OpenMCâs example Pincell Depletion.
Ok, letâs begin
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 equallethargy energies to put in filter
energies = np.logspace(np.log10(1e5), 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 = ['nufission'] # 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() # [neutronscm/source]
flux_unc = t.std_dev.ravel()
# The unit of flux is [n/cmsrc], we need to convert it into [n/cm2s]
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.2e11 # [J/fission]
nu = neu_fission/fiss_rate # [neutrons/fission]
k = neu_fission # [neutrons/source]
Flux = flux/V*(P*nu/(Q*k)) # [neutrons/cm^2sec]
# 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 [neutroncm/source], we need to convert it to [neutron/cm^2sec]

flux â [neutroncm/source]

flux/V = [neutron/cm^2source]
P â reactor power of your element [J/s]
Q = 200 MeV for U235 or 3.2x1011 [J/fission]
nu â score ânufissionâ/'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^2sec]

Also, I gave the code about how to get the Normalized flux/lethargy width.
Some detailed discussions or documents about this topic.
 OpenMCâs Document about Normalization of tally results
 Question about power designation while depleting

OpenMCâs Normalization Factor (This one gives an detailed example)
 Normalization Factor in Section 8.3
 Unit transformation in a fixed source problem
 To find neutron flux, atomic density, and burnup 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.