Developing Axial Flux Distribution (Fast, Epi, Thermal) on a hexagonal Core

I am trying to simulate the design of the core with a Hexagonal Geometry. Can you please describe to me how to develop the axial distribution of flux in thermal, epithermal and fast regions. I am new to openmc and can’t get myself familiar with the type of meshes or tallies to be used, especially in this case. Also for the radial thermal flux distribution. Currently it’s giving me an error that maximum number of lost particles have been reached. I have tried using the spatial legendre filter, please let me know if my approach is correct or not.

Also please if you can guide me on how to learn to use different settings and tallies, since examples help with understanding, not just theoretical description, for things like how to develop plots for peaking factors against radial or axial directions, axial and thermal distributions of flux or power density distributions, especially ranging specific locations like at hot plate, at bottom plate, etc. all using the python api.

Hello! The error you are getting (maximum number of lost particles has been reached) suggests that there’s a problem with your geometry. Check your geometry’s boundary conditions and cell definitions to fix whatever the issue might be that is causing particles to leak out.

As for the axial flux distribution, define a rectilinear mesh with fixed cell sizes from the bottom to the top of your geometry. Then, define a tally that uses that mesh, with an energy filter that has 3 bins: one for thermal, epithermal, and fast neutrons. It should look like this:

mesh = openmc.RegularMesh()
mesh.dimension = [1, 1, 30] # for 30 axial cells
mesh.lower_left = [-0.1, 0.1, -150] # depending on the dimensions of your fuel
mesh.width = [0.1, 0.1, 10] # depending on the dimensions of your fuel

mesh_filter = openmc.MeshFilter(mesh)
energy_filter = openmc.EnergyFilter([0.0, 4.0, 100, 1.0e6])
t = openmc.Tally(tally_id=1)
t.filters.append(mesh_filter, energy_filter)
t.scores = [‘flux’]
tallies = openmc.Tallies([t])
tallies.export_to_xml()

1 Like

Hi there,
Thank you for the answer.

Can you please tell me how I can get this flux value to be normalized against 100MWth power of my reactor?

Can you please guide me on a few things here, my fuel pin pitch is 1.38 cm while assembly pitch is 16 and core pitch is found to be 145.08. I am getting the following errors (as seen in the pics attached). I don’t happen to be an avid python user, can you please correct me where I am going wrong.
Also this core that is designed in 2D and only in infinite z dimension, to apply the mesh I have used the core height of 285cm. But using the tally I am getting errors (also attached for reference).
Screenshot (143)
Do these errors correspond to some mistake in my geometry? Can you help me identify if any, it would be a great help.

Screenshot (144)
This is the respective error, I am getting.

@Issacc The filters attribute of a tally object is a list of filters, so the line where your code is failing should either be:

t.filters.extend([mesh_filter, energy_filter])

or

t.filters = [mesh_filter, energy_filter]
2 Likes

@Issacc You can use f_R = (P_r)_max/ average of P_R formula for radial PPF.

  • code
tally_fn = openmc.Tally(name= 'radial peaking')
tally_fn.filters = [openmc.DistribcellFilter([fuel_cell])]
tally_fn.scores = ['fission']

sp = openmc.StatePoint('statepoint.100.h5')
tally_v = sp.get_tally(name='radial peaking')
fission_rates = tally_v.get_values(scores=['fission'])

fission_rates /= np.mean(fission_rates[fission_rates > 0.])
import numpy as np
np.max(fission_rates)

For axial PPF, f_Z = (P_z)_max/average of P_Z.

  • code
mesh_axial = openmc.RegularMesh(mesh_id=1)
mesh_axial.dimension = [1, 1, 10]
mesh_axial.lower_left = [x1, y1, z1]
mesh_axial.upper_right = [x2, y2, z2]
mesh_f_axial = openmc.MeshFilter(mesh_axial)

tally_axial = openmc.Tally(name='axial peaking')
tally_axial.filters = [mesh_f_axial]
tally_axial.scores = ['fission']

sp = openmc.StatePoint('statepoint.100.h5')
tally_v = sp.get_tally(name='axial peaking')
fission_rates = tally_v.get_values(scores=['fission'])

fission_rates /= np.mean(fission_rates[fission_rates > 0.])
np.max(fission_rates)
1 Like