Trying to model SD-TMSR using OpenMC

I am new to OpenMC and i am trying to recreate an SD-TMSR model based on this paper here (A neutronic evaluation of a thorium-based molten salt breeder reactor - ScienceDirect).

All the fuel composition, temp., density, geometry characteristics, etc. are the same but when i run the initial eigenvalue calculation, it does not even reach to 1 (only ~0.5). The leakage fraction is low (0.000xx), the absorption tally on fuel filters appear to be 99.2%.

Any help? Thanks:D

Could you share your input file?

The low keff value are solved yesterday because I forgot to enrich the Li7 to 99.995% and so Li6 absorbed almost all the neutrons.

But I have another 2 problems that the keff differs when I use different material definition methods for LiF-BeF2-ThF4-UF4 with 70-17.5-12.3-0.2 mole%.

fuel = openmc.Material(name=‘LiF-BeF2-ThF4-UF4’)
lif = openmc.Material(name=‘LiF’)
lif.add_nuclide(‘Li7’, 0.99995, ‘ao’)
lif.add_nuclide(‘Li6’, 0.00005, ‘ao’)
lif.add_element(‘F’, 1, ‘ao’)
bef2 = openmc.Material(name=‘BeF2’)
bef2.add_element(‘Be’, 1)
bef2.add_element(‘F’, 2)
uf4 = openmc.Material(name=‘UF4’)
uf4.add_nuclide(‘U233’, 1)
uf4.add_element(‘F’, 4)
thf4 = openmc.Material(name=‘ThF4’)
thf4.add_nuclide(‘Th232’, 1)
thf4.add_element(‘F’, 4)
fuel = openmc.Material.mix_materials([lif, bef2, thf4, uf4], [0.280403, 0.126674, 0.583406, 0.009517], ‘wo’)
got keff = 0.998 in pin-only sim

fuel = openmc.Material(name=‘LiF-BeF2-ThF4-UF4’)
fuel.add_nuclide(‘Li6’, 3.2417710**(-6), ‘wo’)
fuel.add_nuclide(‘Li7’, 7.5619
10**(-2), ‘wo’)
fuel.add_nuclide(‘F19’, 2.047810**(-1), ‘wo’)
fuel.add_nuclide(‘Be9’, 2.4285
10**(-2), ‘wo’)
fuel.add_nuclide(‘F19’, 1.023910**(-1), ‘wo’)
fuel.add_nuclide(‘U233’, 7.1768
10**(-3), ‘wo’)
fuel.add_nuclide(‘F19’, 2.340310**(-3), ‘wo’)
fuel.add_nuclide(‘Th233’, 4.3948
10**(-1), ‘wo’)
fuel.add_nuclide(‘F19’, 1.4393*10**(-1), ‘wo’)

The weight fraction reference is based on the reference and the calculation I did. got keff - 0.15 on pin-only sim.

Another problem is that while using the first material definition method, I got keff = 0.98 for the full core, and 1.01 for the inner pin only. I found out that when I simulate only the outer pin, the result showed subcriticality (0.82) that I assumed was caused by the fuel-to-moderator ratio that isn’t ideal, so when I changed the core only to 1 region, the criticality was achieved.

However since I am doing benchmarking simulation from the paper, I shouldn’t have to modify the geometry or the material composition.

Here’s the geometry code:

Inner pin fuel channel

ir_pin = openmc.ZCylinder(r=3.5)
ir_mode = openmc.model.HexagonalPrism(edge_length=7.5, boundary_type = ‘reflective’)
itop_core = openmc.ZPlane(z0=230.0, boundary_type = ‘reflective’)
ibot_core = openmc.ZPlane(z0=-230.0,boundary_type = ‘reflective’)
ifuel_cell = openmc.Cell(fill=fuel, region=-ir_pin & -itop_core & +ibot_core)
imoderator_cell = openmc.Cell(fill=moderator, region=+ir_pin & -ir_mode & -itop_core & +ibot_core)
ipin_universe = openmc.Universe(cells=[ifuel_cell, imoderator_cell,])

Outer pin fuel channel

or_pin = openmc.ZCylinder(r=5.0)
or_mode = openmc.model.HexagonalPrism(edge_length=7.5, boundary_type = ‘reflective’)
otop_core = openmc.ZPlane(z0=215.0, boundary_type = ‘reflective’)
obot_core = openmc.ZPlane(z0=-215.0,boundary_type = ‘reflective’)
ofuel_cell = openmc.Cell(fill=fuel, region=-or_pin & -otop_core & +obot_core)
omoderator_cell = openmc.Cell(fill=moderator, region=+or_pin & -or_mode & -otop_core & +obot_core) #
opin_universe = openmc.Universe(cells=[ofuel_cell, omoderator_cell,])

geometry = openmc.Geometry(opin_universe)
geometry.export_to_xml()

batches = 100
inactive = 30
particles = 2000
settings_file = openmc.Settings()
settings_file.batches = batches
settings_file.inactive = inactive
settings_file.particles = particles

openmc.run()

So, I have a couple suggestions from what I see here. I think one issue may be the number of particles being low, so you might have larger than expected uncertainty in your keff value. Another note is that the density isn’t set in the material definitions you use, so if you check the materials.xml fuel density it probably isn’t the density you want. If you add fuel.set_density('g/cm3', 3.3) then that should fix that (if you’re using 3.3 g/cc for your density). For the second material definition, I’m not certain if OpenMC reads in ‘F19’ correctly if it’s defined multiple times with different weight percentages (though it might), but you could consider changing that to a single line of fuel.add_nuclide('F19', 4.5e-1) to see if that changes anything. Hope that helps!