OpenMC Dose Calculation: How to Determine Optimal Settings for Accuracy

User Support
1

Hello everyone,

I am currently running an OpenMC simulation focused on neutron flux and dose calculations from a fixed neutron source (3e10 n/sec) located in the center of a room with concrete walls (2D geometry). I have a few questions regarding the simulation setup and would greatly appreciate your help.

  1. How can I determine the sensitivity of the simulation results to the number of particles and the number of batches I set? What factors should I consider when selecting these parameters?

  2. Should I assess the effect of mesh resolution on the results? If so, how can I check its influence on the accuracy of the flux and dose calculations?

  3. Regarding the flux tally, if I understand correctly, the calculation is normalized by the source particle flux (particles/sec). So, in my case, to calculate the neutron flux and dose for my specific problem, I would need to multiply the result by the neutron source rate (3e10 n/sec). Is this the correct approach?

  4. Finally, assuming the above is correct, how should I verify the effect of the neutron source strength on the dose calculations? Are there any recommended methods for this verification?

Here is a part of my code that I am working on:

"source_position = (xplane_2.x0 + width_inside_room/2, yplane_2.y0 + depth_inside_room/2, 0)
space = openmc.stats.Point(source_position)
angle = openmc.stats.Isotropic()
energy = openmc.stats.Discrete([14e6], [1.0])
source = openmc.IndependentSource(space=space, angle=angle, energy=energy)
source.particle = “neutron”

# Mesh and filter definitions
mesh = openmc.RegularMesh().from_domain(geometry)
mesh.dimension = [800, 800]
mesh_filter = openmc.MeshFilter(mesh)
particle_filter = openmc.ParticleFilter('neutron')

# Flux tally
flux_tally = openmc.Tally(name="flux tally")
flux_tally.filters = [mesh_filter, particle_filter]
flux_tally.scores = ["flux"]

# Dose tally
energy_bins_n, dose_coeffs_n = openmc.data.dose_coefficients(particle="neutron", geometry="AP")
energy_function_filter_n = openmc.EnergyFunctionFilter(energy_bins_n, dose_coeffs_n)
energy_function_filter_n.interpolation = "cubic"

dose_tally = openmc.Tally(name="neutron_dose_on_mesh")
dose_tally.filters = [mesh_filter, particle_filter, energy_function_filter_n]
dose_tally.scores = ["flux"]

tallies = openmc.Tallies([flux_tally, dose_tally])

# Settings
settings = openmc.Settings()
settings.run_mode = "fixed source"
settings.source = source
settings.particles = 10000000
settings.batches = 100

# Create and run model
model = openmc.Model(geometry, materials, settings, tallies)
sp_filename = model.run()

# Process results
 flux_tally = sp.get_tally(name="flux tally")
  dose_tally = sp.get_tally(name="neutron_dose_on_mesh")

    flux_mean = flux_tally.get_reshaped_data(value='mean', expand_dims=True).squeeze()
    flux_mean = flux_mean * neutron_flux
    flux_std_dev = flux_tally.get_reshaped_data(value='std_dev', expand_dims=True).squeeze()
    dose_mean = dose_tally.get_slice(scores=['flux']).mean.reshape(mesh.dimension) "

Thank you very much for your guidance.

2

For this we know the answer without variance reduction the uncertainty on the tally will be proportional to \frac{1}{\sqrt{n}}, where n is the number of histories ran. In this case you will be using a fixed source problem, and so batches don’t really matter (that’s only an eigenvalue thing) and just the total number of particles ran matters. You can do a bunch of back of the envelope, or you can extrapolate. Just run a small number of histories on the order of 10k - 100k. Look at the tally uncertainties, and then extrapolate to get the desired uncertainty. If you have 20% uncertainty and want 2%, then run 100x number of histories.

Yes, in general this is a good idea to do. Though for such a simple geometry you probably won’t a super fine mesh. You will see tradeoff between mesh volume and uncertainty. The smaller the mesh the more histories you need to get the same uncertainties. You should try to use the coarsest mesh for resolving the phenomenon you are trying to model.

Yes, you are correct. This is one of the easiest problems to calculate the number of source particles, so you lucked out there. Just remember that your dose rate will be in seconds, and not hours.

I’m a fan of hand calculations. I should have the equations to do this by hand off the top of my head, but I don’t. If your dose rate is within 20% or so there probably wasn’t an error.

3

Thank you so much for your help!

I have another question and was wondering if you could please help me with it. Would the leakage fraction result output provide useful information about the simulation and its accuracy?"

4

My understanding of the leaking fraction is that in this case: not really. All leakage will tell you is how many particles are making it to the vacuum boundary of your problem, which is completely reasonable to happen in this case. If you had reflective boundaries the leakage should be 0.

5

Thank you so much! I really appreciate your help.

I have another question regarding photon dose calculation and was wondering if you could help me with it. I am trying to calculate the dose from a Co-60 point source, which is shielded by lead boundaries, and further by concrete walls. Here is the code I have for the dose calculation part:

Defining Co60 source

strength_per_rod = 7500 * 3.7e10 # 7500 Ci converted to Bq
total_strength = 3 * strength_per_rod # 3 rods
source = openmc.Source()
source.space = openmc.stats.Point((outside_width + concrete_wall_width + Co60_distance_left_wall, outside_depth + concrete_wall_depth + floor_plan_depth/2, height/2))
source.angle = openmc.stats.Isotropic()
source.energy = openmc.stats.Discrete([1.17e6, 1.33e6], [0.5, 0.5])
source.particle = “photon”
source.strength = total_strength

Settings

settings = openmc.Settings()
settings.output = {“tallies”: False}
settings.batches = 80
settings.inactive = 10
settings.particles = 100000000
settings.run_mode = ‘fixed source’
settings.source = source
settings.photon_transport = True

Create mesh tally for photon flux

mesh = openmc.RegularMesh()
mesh.dimension = [100, 100, 1]
mesh.lower_left = [xplane_0.x0, yplane_0.y0, zplane_0.z0]
mesh.upper_right = [xplane_12.x0, yplane_11.y0, zplane_1.z0]
mesh_filter = openmc.MeshFilter(mesh)
total_width = xplane_12.x0 - xplane_0.x0
total_depth = yplane_11.y0 - yplane_0.y0
total_height = zplane_1.z0 - zplane_0.z0
mesh_element_volume = (total_width / mesh.dimension[0]) * (total_depth / mesh.dimension[1]) * (total_height / mesh.dimension[2])
particle_filter = openmc.ParticleFilter(‘photon’)

Create mesh tally for calculating Photon flux

flux_tally = openmc.Tally(name=‘flux tally’)
flux_tally.filters = [mesh_filter, particle_filter]
flux_tally.scores = [‘flux’]

Dose tally

energy_bins_n, dose_coeffs_n = openmc.data.dose_coefficients(particle=“photon”, geometry=“AP”)
energy_function_filter_n = openmc.EnergyFunctionFilter(energy_bins_n, dose_coeffs_n)
energy_function_filter_n.interpolation = “cubic”
dose_tally = openmc.Tally(name=“photon_dose_on_mesh”)
dose_tally.filters = [mesh_filter, particle_filter, energy_function_filter_n]
dose_tally.scores = [“flux”]
tallies = openmc.Tallies([flux_tally, dose_tally])
model = openmc.model.Model(geometry, materials, settings, tallies)
sp_filename = model.run()
with openmc.StatePoint(sp_filename) as sp:
flux_tally = sp.get_tally(name=“flux tally”) # flux unit in OpenMC is: Particle-cm/source particle
dose_tally = sp.get_tally(name=“photon_dose_on_mesh”) # dose unit in OpenMC is: Psv-cm^3
flux_mean = flux_tally.get_reshaped_data(value=‘mean’, expand_dims=True).squeeze()
flux_std_dev = flux_tally.get_reshaped_data(value=‘std_dev’, expand_dims=True).squeeze()
dose_mean = dose_tally.get_slice(scores=[‘flux’]).mean.reshape(mesh.dimension)
dose_std_dev = dose_tally.get_slice(scores=[‘flux’]).std_dev.reshape(mesh.dimension)

convert flux to particle/cm^2/s

flux_mean = (flux_mean/mesh_element_volume) * total_strength
flux_std_dev = (flux_std_dev/mesh_element_volume) * total_strength

Calculate relative error for flux

flux_relative_error = np.divide(flux_std_dev, flux_mean, out=np.zeros_like(flux_std_dev), where=flux_mean!=0)

convert dose to rem/year

total_dose = (dose_mean/mesh_element_volume) * total_strength # Psv/sec
total_dose = total_dose * 1e-9 # mSv/sec
total_dose = total_dose * 60 * 60 * 24 * 364 # mSv/year
total_dose = total_dose * 0.1 # rem/year

Convert dose standard deviation to rem/year

dose_std_dev = (dose_std_dev/mesh_element_volume) * total_strength # Psv/sec
dose_std_dev = dose_std_dev * 1e-9 # mSv/sec
dose_std_dev = dose_std_dev * 60 * 60 * 24 * 364 # mSv/year
dose_std_dev = dose_std_dev * 0.1 # rem/year

Calculate relative error for dose

dose_relative_error = np.divide(dose_std_dev, total_dose, out=np.zeros_like(dose_std_dev), where=total_dose!=0)

Print some statistics

print(f"Max flux relative error: {np.max(flux_relative_error):.2%}“)
print(f"Mean flux relative error: {np.mean(flux_relative_error):.2%}”)
print(f"Max dose relative error: {np.max(dose_relative_error):.2%}“)
print(f"Mean dose relative error: {np.mean(dose_relative_error):.2%}”)

However, when I calculate the standard deviation and relative error of the results, I get very high values, especially near the point source and up to the shielding walls. I tried adjusting the mesh resolution and increasing the number of particles, but this did not improve the results. I have not used any variance reduction method. Do you have any suggestions for me? Alternatively, is there a more precise method for calculating the photon dose that you would recommend?

Many thanks!