Dose calculation around a punctual source

Hello !

I have a misunderstanding with the dose calculation on the surface of a sphere around a point source.

I made several tests by varying the density of the material of the sphere containing the punctual source (‘sic_core’).
I expected the dose to decrease with increasing density (self-shielding).
On the contrary, the dose increased. Therefore, I wonder if I am measuring a dose…

My code is the following:

import openmc
import math


# Define Materials
# Air
air = openmc.Material(name='Air')
air.set_density('g/cm3', 0.001205)
air.add_element('N', 0.784431)
air.add_element('O', 0.210748)
air.add_element('Ar',0.0046)
# SiC with reduced density
sic_core = openmc.Material(name='SiC_core')
sic_core.set_density('g/cm3', 3.15*0.22)
sic_core.add_element('Si', 0.5, 'ao')
sic_core.add_element('C', 0.5, 'ao')

materials = openmc.Materials([air, sic_core])


# Define Geometry
# rayons
radius_calcul_cm = 1
radius_exter_cm = 2
# surfaces
surface_calcul = openmc.Sphere(r=radius_calcul_cm)
surface_exter = openmc.Sphere(r=radius_exter_cm, boundary_type='vacuum')
# cells
interior_cell = openmc.Cell(fill = sic_core, region= -surface_calcul)
exterior_cell = openmc.Cell(fill = air, region= +surface_calcul & -surface_exter)
# universe
universe = openmc.Universe(cells=[interior_cell, exterior_cell])
# geometry
simple_geometry = openmc.Geometry(universe)


# Define Sources
source = openmc.Source()
source.space = openmc.stats.Point()
source.angle = openmc.stats.Isotropic()
source.energy = openmc.stats.Discrete([1e+6], [1.0])


# Define Settings
sett = openmc.Settings()
sett.photon_transport = True
sett.batches = 100
sett.inactive = 0
sett.particles = 10000
sett.run_mode = 'fixed source'
sett.source = source


# Define Tally
# Filters
# particule filter
photon_particle_filter = openmc.ParticleFilter(['photon'])
# dose filter
energies, dose_coeffs = openmc.data.dose_coefficients('photon', 'AP')
print(energies, dose_coeffs)
dose_filter = openmc.EnergyFunctionFilter(energies, dose_coeffs)
# surface filter
surface_calcul_filter = openmc.SurfaceFilter(surface_calcul)
# Tallies
tallies = openmc.Tallies()
# Dose in contact with the SiC
dose_contact_tally = openmc.Tally(name='dose_contact')
dose_contact_tally.filters = [
    photon_particle_filter,
    surface_calcul_filter,
    dose_filter
]
dose_contact_tally.scores = ['current']

tallies.append(dose_contact_tally)

# Define Model
model = openmc.model.Model(
    geometry=simple_geometry,
    materials= materials,
    settings= sett,
    tallies= tallies)

results_filename = model.run()

Hi @Arthurcrrd! Welcome to the community!

I haven’t had a chance to run this model myself, but from a quick glance at the geometry, materials, and source I’m guessing that most if not all of the particles are reaching the tally surface with few to no interactions for the nominal density of silicon carbide. Because this is a neutron source and a “photon” particle filter is applied to the tally, there are going to be very few secondary photons produced from neutron reactions, resulting in a lower dose. As the density of the material increases, you’re right that there will be more self-shielding of the neutrons but the number of photons generated will increase and in turn the tallied photon dose will increase.

As an exercise, if you change the source term from a neutron source to a photon source you’ll see the dose rate change with the behavior you’ve described above.

Hi @pshriwise !
Thank you for the answer !

Yes, you’re right.
I just realized that my source was not necessarily photons but neutrons by default!
My intention was to put a photon fixed source, not neutron.
The results I have are logical. I’m just calculating the secondary gamma dose.

So I have to add the line : source.particule = 'photon'