Dear OpenMC users,
I am trying to calculate the thermalization efficiency for the HDPE moderator which is 7.5 cm thick. I am placing two small silver foils (1 cm by 1 cm by 0.02 cm) in front and behind the HDPE block to record flux values on them. Below, I have attached the main part of my code. I am using a flux tally with an energy filter to calculate the thermalized neutron flux after the moderator. The efficiency I obtain is around 1%, which is lower than what I have seen in the literature. I was wondering if there is anything I might be overlooking in my code, or if I should use a different approach for this simulation. If you know of any examples that could help, it would be much appreciated. Many thanks for your help.
Material definitions
mat_air = openmc.Material(name=“Air”)
mat_air.add_element(“N”, 0.784431)
mat_air.add_element(“O”, 0.210748)
mat_air.add_element(“Ar”, 0.0046)
mat_air.set_density(“g/cc”, 0.001205)
mat_hdpe = openmc.Material(name = “HDPE”)
mat_hdpe.add_element(“C”, 2)
mat_hdpe.add_element(“H”, 4)
mat_hdpe.set_density(“g/cm3”, 0.95)
mat_hdpe.add_s_alpha_beta(‘c_H_in_CH2’)
aw_Ag = 107.87 # silver
aw_Cu = 63.55 # Copper
mass_frac_Ag = 0.925
mass_frac_Cu = 0.075
mole_Ag = mass_frac_Ag/aw_Ag
mole_Cu = mass_frac_Cu/aw_Cu
total_moles = mole_Ag + mole_Cu
atom_frac_Ag = mole_Ag/total_moles
atom_frac_Cu = mole_Cu/total_moles
mat_sterling_silver = openmc.Material(name = “Sterling Silver”)
mat_sterling_silver.add_element(‘Ag’, atom_frac_Ag)
mat_sterling_silver.add_element(‘Cu’, atom_frac_Cu)
mat_sterling_silver.set_density(“g/cm3”, 9.70)
materials = openmc.Materials([mat_air, mat_hdpe, mat_sterling_silver])
neutron_flux = 1e9
space = openmc.stats.Point((x_source, y_source, z_source))
angle = openmc.stats.Isotropic()
energy = openmc.stats.Discrete([2.45e6], [1])
source = openmc.IndependentSource(space = space, angle = angle, energy = energy)
source.particle = “neutron”
Total flux at foil #1
particle_filter = openmc.ParticleFilter(‘neutron’)
total_flux_tally_foil1 = openmc.Tally(name = “total flux tally foil 1”)
total_flux_tally_foil1.filters = [openmc.CellFilter(silver_foil1_cell), particle_filter]
total_flux_tally_foil1.scores = [“flux”]
Total flux at foil #2
total_flux_tally_foil2 = openmc.Tally(name = “total flux tally foil 2”)
total_flux_tally_foil2.filters = [openmc.CellFilter(silver_foil2_cell), particle_filter]
total_flux_tally_foil2.scores = [“flux”]
Thermalized neutron flux at foil #2
thermalized_flux_tally_foil2 = openmc.Tally(name = “thermalized flux tally foil 2”)
energy_filter = openmc.EnergyFilter([0.0, 0.1]) # Energy group below 0.05 eV
thermalized_flux_tally_foil2.filters = [openmc.CellFilter(silver_foil2_cell), particle_filter, energy_filter]
thermalized_flux_tally_foil2.scores = [“flux”]
tallies = openmc.Tallies([total_flux_tally_foil1, total_flux_tally_foil2, thermalized_flux_tally_foil2])
Define setting
settings = openmc.Settings()
settings.output = {“tallies”: False}
settings.run_mode = “fixed source”
settings.source = source
settings.batches = 35
settings.inactive = 0
settings.particles = 100000000
settings.geometry_debug = True
Process results
with openmc.StatePoint(sp_filename) as sp:
total_flux_tally_foil1_result = sp.get_tally(name = “total flux tally foil 1”) # flux unit in OpenMC is particle-cm/source
total_flux_tally_foil1_mean = total_flux_tally_foil1_result.get_reshaped_data(value=‘mean’, expand_dims=True).squeeze()
total_flux_tally_foil1_std_dev = total_flux_tally_foil1_result.get_reshaped_data(value=‘std_dev’, expand_dims=True).squeeze()
total_flux_tally_foil1_mean = total_flux_tally_foil1_mean/Silver_foil_volume
total_flux_tally_foil1_std_dev = total_flux_tally_foil1_std_dev/Silver_foil_volume
total_flux_tally_foil2_result = sp.get_tally(name = “total flux tally foil 2”) # flux unit in OpenMC is particle-cm/source
total_flux_tally_foil2_mean = total_flux_tally_foil2_result.get_reshaped_data(value=‘mean’, expand_dims=True).squeeze()
total_flux_tally_foil2_std_dev = total_flux_tally_foil2_result.get_reshaped_data(value=‘std_dev’, expand_dims=True).squeeze()
total_flux_tally_foil2_mean = total_flux_tally_foil2_mean/Silver_foil_volume
total_flux_tally_foil2_std_dev = total_flux_tally_foil2_std_dev/Silver_foil_volume
thermalized_flux_tally_foil2_result = sp.get_tally(name = “thermalized flux tally foil 2”) # flux unit in OpenMC is particle-cm/source
thermalized_flux_tally_foil2_mean = thermalized_flux_tally_foil2_result.get_reshaped_data(value=‘mean’, expand_dims=True).squeeze()
thermalized_flux_tally_foil2_std_dev = thermalized_flux_tally_foil2_result.get_reshaped_data(value=‘std_dev’, expand_dims=True).squeeze()
thermalized_flux_tally_foil2_mean = thermalized_flux_tally_foil2_mean/Silver_foil_volume
thermalized_flux_tally_foil2_std_dev = thermalized_flux_tally_foil2_std_dev/Silver_foil_volume