Conventional PWR Pincell Code giving a fast flux

Hi there,

I needed to compare my reactor’s neutron spectrum with a conventional pwr pincell so I used the pincell program in examples present as the legacy codes when installing openmc. I required values for the scoping analysis and the neutron ratios in thermal, epithermal and fast regimes showed me a fast spectrum. Is it that I am having a bug or am I making any mistakes in doing so.

Following is the code
Material specifications:
(I haven’t added any boron reactivity, since I need the comparison without it, and the material specs are the same as I need for my case)

fuel = openmc.Material(name='Fuel')

helium = openmc.Material(name='Helium for gap')
helium.set_density('g/cm3', 0.001598)
helium.add_element('He', 2.4044e-4)

clad = openmc.Material(name='Clad')
clad.add_nuclide('Cr50', 3.30121e-6)
clad.add_nuclide('Cr52', 6.36606e-5)
clad.add_nuclide('Cr53', 7.2186e-6)
clad.add_nuclide('Cr54', 1.79686e-6)
clad.add_nuclide('Fe54', 8.68307e-6)
clad.add_nuclide('Fe56', 1.36306e-4)
clad.add_nuclide('Fe57', 3.14789e-6)
clad.add_nuclide('Fe58', 4.18926e-7)
clad.add_nuclide('Zr90', 2.18865e-2)
clad.add_nuclide('Zr91', 4.77292e-3)
clad.add_nuclide('Zr92', 7.29551e-3)
clad.add_nuclide('Zr94', 7.39335e-3)
clad.add_nuclide('Zr96', 1.1911e-3)
clad.add_nuclide('Sn112', 4.68066e-6)
clad.add_nuclide('Sn114', 3.18478e-6)
clad.add_nuclide('Sn115', 1.64064e-6)
clad.add_nuclide('Sn116', 7.01616e-5)
clad.add_nuclide('Sn117', 3.70592e-5)
clad.add_nuclide('Sn118', 1.16872e-4)
clad.add_nuclide('Sn119', 4.14504e-5)
clad.add_nuclide('Sn120', 1.57212e-4)
clad.add_nuclide('Sn122', 2.23417e-5)
clad.add_nuclide('Sn124', 2.79392e-5)
clad.add_nuclide('Hf174', 3.54138e-9)
clad.add_nuclide('Hf176', 1.16423e-7)
clad.add_nuclide('Hf177', 4.11686e-7)
clad.add_nuclide('Hf178', 6.03806e-7)
clad.add_nuclide('Hf179', 3.0146e-7)
clad.add_nuclide('Hf180', 7.76449e-7)
clad.temperature = 600

coolant = openmc.Material(name='Coolant')
coolant.set_density('g/cm3',, pressure=15.5132))
coolant.add_element('H', 2.0)
coolant.add_element('O', 1.0)
coolant.temperature = 600

# Collect the materials together and export to XML
materials = openmc.Materials([fuel, helium, clad, coolant])

Geometry specifications:

# Create cylindrical surfaces
fuel_or = openmc.ZCylinder(r=0.39218, name='Fuel OR')
clad_ir = openmc.ZCylinder(r=0.40005, name='Clad IR')
clad_or = openmc.ZCylinder(r=0.45720, name='Clad OR')

# Create a region represented as the inside of a rectangular prism
pitch = 1.25984
box = openmc.rectangular_prism(pitch, pitch, boundary_type='reflective')

# Create cells, mapping materials to regions
fuel = openmc.Cell(fill=fuel, region=-fuel_or)
gap = openmc.Cell(fill=helium, region=+fuel_or & -clad_ir)
clad = openmc.Cell(fill=clad, region=+clad_ir & -clad_or)
water = openmc.Cell(fill=coolant, region=+clad_or & box)

# Create a geometry and export to XML
geometry = openmc.Geometry([fuel, gap, clad, water])


settings = openmc.Settings()
settings.run_mode= "eigenvalue"
settings.batches = 120
settings.inactive = 10
settings.particles = 5000

And the tallies to get the flux values:

energies = [0.0, 0.625, 5.531e3, 20e6]
energy_filter = openmc.EnergyFilter(energies)
spectrum_tally = openmc.Tally(name="Flux spectrum")
spectrum_tally.filters = [energy_filter]
spectrum_tally.scores = ['flux']

mesh = openmc.RegularMesh()
mesh.dimension = (1, 1)
mesh.lower_left = (-pitch/2.0, -pitch/2.0)
mesh.upper_right = (pitch/2.0, pitch/2.0)
# Create a mesh filter that can be used in a tally
mesh_filter = openmc.MeshFilter(mesh)

mesh_tally = openmc.Tally(name="Mesh tally")
mesh_tally.filters = [mesh_filter]
mesh_tally.scores = ['flux', 'fission', 'nu-fission']

tallies = openmc.Tallies([spectrum_tally, mesh_tally])

And this is what it plots:
Screenshot (230)

And these are the respective results I am getting (to normalize the flux I have used the power specifications from beavrs for a 3411MWth reactor):

I would really appreciate if someone might point out where I am doing what wrong, because this gives me a hard spectrum for a pwr pincell and the fractions are also not satisfactory.


There is actually nothing wrong in your results. Neutrons are born from fission at fast energies and eventually scatter down to thermal energies. It is expected that in an LWR you will have an appreciable flux in the epithermal and fast energies. For example, see the tallied flux spectrum at the end of this example notebook.

1 Like

Thank you for the response.
Can you please tell me if this is the nature of flux profile I will get even for a PWR, how can I then find out whether a reactor under analysis, is a fast reactor or a thermal one. Is there any other analysis that I can do to find out all the same? I was comparing it with this PWR to figure this out.

It might be harder to reason about the spectrum when you only have a a few energy groups. If you plot the flux with hundreds of energy groups, it becomes a little more obvious what you’re looking at because you would see a typical fission spectrum at high energies, a Maxwellian distribution at low energies where neutrons have thermalized, and then a slowing-down spectrum in between. In a fast reactor, you won’t see the Maxwellian at low energies at all.

Perhaps a better approach here would be to look at the average energy causing fission – this could be done by tallying the fission reaction rate with an energy filter and then finding the average of the energy multiplied by fission rate.

Just to add there are some nice energy GROUP_STRUCTURES built in to OpenMC

Here is an example using VITAMIN-J-175 but there are other options

#creates an empty tally object
tallies = openmc.Tallies()

# sets up filters for the tallies
neutron_particle_filter = openmc.ParticleFilter(['neutron'])
energy_bins = openmc.mgxs.GROUP_STRUCTURES['VITAMIN-J-175']
energy_filter = openmc.EnergyFilter(energy_bins)

# setup the filters for the cell tally
cell_filter = openmc.CellFilter(cell_name_here) 

# create the tally
cell_spectra_tally = openmc.Tally(name='cell_spectra_tally')
cell_spectra_tally.scores = ['flux']
cell_spectra_tally.filters = [cell_filter, neutron_particle_filter, energy_filter]