Different k-eff in RESULTS and calculation

Hi, everyone!

I’ve tried to calculate the material arrangement of three concentric circles on a two-dimensional plane, which are fuel, monolith and reflector respectively from inside to outside. The boundary condition is the white boundary. And I use 50,0000 particles for 300 batches (100 inactive batches).

However, the k-eff shown in RESULTS is 1.86732 +/- 0.00009, and the k-eff calculated from fiss_rate / abs_rate (since the leak is zero) is 1.720066 +/- 0.000238.

I don’t know how this happens. Please excuse my lack of knowledge on OpenMC as I am practically new to it and still on a learning curve.

Thank you for your time and effort!

@Mio welcome to the community! k-effective is defined as production / loss, and the fission rate alone doesn’t tell you the production; for that you need the “nu-fission” tally score, which accounts for the number of neutrons produced per fission (nu). Once you account for that, you should get a number closer to what is reported for k-effective.

Thank you for your help! But actually I did use the “nu-fission” tally score in my test.

Here are some of the codes for K-Eigenvalue (infinity) tallies:

tallies_file = openmc.Tallies()
fiss_rate = openmc.Tally(name=‘fiss. rate’)
abs_rate = openmc.Tally(name=‘abs. rate’)
fiss_rate.scores = [‘nu-fission’]
abs_rate.scores = [‘absorption’]
tallies_file += (fiss_rate, abs_rate)

Hi, Paul. I still can’t know what’s going wrong in my code. Could you please help me?
Since I’m a new user and I can’t upload my .py file, I’ll paste all the codes below.

# import
import openmc
# geometry info
Fuel_D = 1.5
Ref_inner = 2.5
Ref_outer = 5.5
particles = 500000
inactive = 100
batches = 300
group_edges = [0., 0.625, 9.11882e3,4.9787e5,20.0e6]
# material settings
fuel = openmc.Material(name='UC')
fuel.temperature= 1000

Basis = openmc.Material(name='Graphite')
Basis.temperature= 1000

Reflector = openmc.Material(name='Be')
Reflector.temperature= 1000

materials_list = [fuel, Basis, Reflector]
materials = openmc.Materials(materials_list[:])
# geometry settings
fuel_wall = openmc.ZCylinder(x0=0 ,y0=0 ,r=Fuel_D/2) #vacuum
Basis_wall = openmc.ZCylinder(x0=0 ,y0=0 ,r=Ref_inner)
Ref_wall = openmc.ZCylinder(x0=0 ,y0=0 ,r=Ref_outer, boundary_type='white') 

root_universe = openmc.Universe(universe_id=123456789, name='root universe')

fuel_cell = openmc.Cell(name='fuel')
fuel_cell.fill = fuel
fuel_cell.region = -fuel_wall

monolith_cell = openmc.Cell(name='monolith')
monolith_cell.fill = Basis
monolith_cell.region = -Basis_wall & +fuel_wall

reflector_cell = openmc.Cell(name='mreflector')
reflector_cell.fill = Reflector
reflector_cell.region = -Ref_wall & +Basis_wall

geometry = openmc.Geometry(root_universe)
# calculation settings
settings = openmc.Settings()
settings.batches = batches
settings.inactive = inactive
settings.particles = particles
settings.temperature['method'] = 'interpolation'

bounds = [-Ref_inner/2,-Ref_inner/2,-Ref_inner/2,Ref_inner/2,Ref_inner/2,Ref_inner/2]
uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:], only_fissionable=True)
settings.source = openmc.source.Source(space=uniform_dist)

# tally settings
tallies_file = openmc.Tallies()

fiss_rate = openmc.Tally(name='fiss. rate')
abs_rate = openmc.Tally(name='abs. rate')
fiss_rate.scores = ['nu-fission']
abs_rate.scores = ['absorption']
tallies_file += (fiss_rate, abs_rate)
# run
# postprocess
sp = openmc.StatePoint('statepoint.{}.h5'.format(batches))

fiss_rate = sp.get_tally(name='fiss. rate')
abs_rate = sp.get_tally(name='abs. rate')

keff = fiss_rate / (abs_rate)


If you look at the result for the absorption score, you’ll see that it is about 1.08. In a problem without leakage, you normally expect to see a value very close to 1 because the units are “per source neutron” and without leakage, each neutron that is born eventually has to be absorbed somewhere. In your problem, you actually have an appreciable amount of (n,2n) reactions occurring, which increases the population. If you add the (n,2n) score to one of your tallies, you’ll see that it is about ~0.08. So, to recreate the k-effective value from these tallies, you can do nu-fission/(absorption - (n,2n)).

There is actually no universal agreement on how to define k-effective when (n,xn) reactions are present. See this 2009 report by Red Cullen that explains how two codes (MCNP and TART) handle it differently. The definition that is used in OpenMC is consistent with MCNP and Serpent and several other codes.

Hi, many thanks for your reply. Your answer has perfectly solved my problem!

But here comes another question. If I want to use OpenMC to generate cross-sections, how to define settings in the tally regarding to this (n,2n) reactions.

Best wishes!

Hi Paul,

Thanks for the answer!

Maybe we should make this clearer: for a problem listed above, we have been trying to use OpenMC to generate multi-group cross sections for deterministic calculations. However, the multi-group cross sections seemed to incur significant error when we compare deterministic results with openMC results. Using the generated multi-group cross-sections and multi-group neutron flux do not reproduce exactly the OpenMC k-effective. As you mentioned, the discrepancy was caused by our omission of (n,xn) reaction rate. In this case, how do we take into account the contribution of (n,xn) reaction when generating multi-group cross-sections?

We tried activating

  1. ‘nu=True’ in mgxs.ScatterXS(), and
  2. ‘nu=True’ in mgxs.TransportXS()
    but it doesn’t help improving the deterministic results. What else should we do to fix this problem?

Thanks for your help.

I’m probably not the best person to address questions about downstream use for deterministic codes (perhaps @mkreher13 or @agnelson can chime in) but I think to get a matching answer, you would need to replace your absorption cross section with a “removal” cross section that accounts for the extra production of neutrons from (n,2n), i.e., removal = absorption - (n,2n). If you want to get an (n,2n) MGXS, you can do so with the ArbitraryXS class.

A few things to keep in mind:

  1. Yes, the removal cross section is a good strategy to address (n,2n) reactions
  2. Use a good group structure for your deterministic method
  3. Beware of diffusion coefficients (if using)

Point 2 is not trivial, depending on your reactor, but there are some good group structures out there for various neutron spectrums.

Point 3 is very important to understand. Monte Carlo tallied diffusion coefficients should be treated with extreme caution.
Please see this discussion: Units of Diffusion_Coefficient? - #5 by bakingbad
and this OpenMC Pull Request to read further: Condensing diffusion instead of transport by mkreher13 · Pull Request #1916 · openmc-dev/openmc · GitHub
Furthermore, if you are using diffusion, consider using a non-linear diffusion coefficient to smooth out the differences between Monte Carlo and your diffusion solver.

Hi Miriam and Paul,

Thank you very much for your timely response. Given your hints, I believe we should try removing the contribution of (n,2n) in the total XS.

By the way, I feel that this could be a underlying trouble for researchers who use OpenMC to generate multi-group XSs. The special treatment of (n,xn) is automatically considered in deterministic lattice codes like Helios and others. And we just found out that this problem was also noticed by SERPENT users:

Revised multi-group XS in SERPENT

Therefore, I thought it might be good to consider adding an homogenization option that accounts for the (n,xn) modification automatically.

Thanks again,


That’s an excellent suggestion @TengfeiZhang. I’ve just created an issue on our GitHub repo to track this.

I’m told that the ‘nu-scatter’ tally accounts for (n,2n) reactions and should be used in this case. See documentation here: 8. Specifying Tallies — OpenMC Documentation

Yes, that’s correct, but I believe we also need to account for it in absorption as well in order for the total cross section to be reproduced accordingly. That is, (n,xn) production is added to scattering and subtracted from absorption so that absorption + scatter = total.

Hi Paul and Miriam, thanks again for the kind help and consideration!

As an update, I’ve just submitted a pull request with a new ReducedAbsorptionXS class.

1 Like