openmc
openmc copied to clipboard
`k_inf` estimator for `IndependentOperator`
It would be nice to be able to estimate $k_\infty$ when using IndependentOperator
so it can can be compared against tools like ORIGEN-S.
There is a method in this paper that seems doable. The estimate is as follows:
$$ k_\infty = \frac{\Sigma_i N_i \Sigma_j \nu_j \sigma_{i,j}}{\Sigma_i N_i \Sigma_j \sigma_{i,j}} \text{ for all } i \text{ nuclides, } j \text{ reactions.} $$
where $N_i$ is the concentration of nuclide $i$, $\sigma_{i,j}$ is the microscopic cross section of reaction $j$ for nuclide $i$, and $\nu_j$ is the neutron multiplication of reaction $j$, defined as
$$ \nu_j = \begin{cases} 0 & j=\text{ absorption reaction}\ x & j=(n,xn)\ \bar{\nu} & j=\text{fission} \end{cases} $$
The most difficult part of getting this formulation would be getting the values for $\bar{\nu}$. One possible path forward would be to figure out a way to energy-average $\nu(E)$ while creating the one-group libraries.
A sample implementation is below. Note that this function is written as if it were declared within IndependentOperator
:
def _estimate_k_inf():
"""Estimate :math:`k_\infty`
:math:`k_\infty = \frac{\Sigma_i N_i \Sigma_j \nu_j \sigma_{i,j}}{\Sigma_i N_i \Sigma_j \sigma_{i,j}}` for all :math:`i` nuclides, :math:`j` reactions.
"""
nucs = self.micro_xs.index
rxns = self.micro_xs.columns
top = 0
bottom = 0
for nuc, rxn in product(nucs, rxns):
# get total nuclide concentration across all materials
ids = np.arange(0, len(self.number.number))
volumes = self.number.get_mat_volume(ids)
N = self.number.get_atom_density(ids, nuc)
N = np.dot(N, volumes) / np.sum(volumes)
sigma = self.micro_xs[rxn].loc[nuc]
re_search = re.search("n,\d*n", rxn)
if bool(re_match):
n = re_search.group(0)[2]
if n == 'n':
n = 1
else:
n = int(n)
nu = n
elif rxn == 'fission':
nu = self._nu_bar
else:
nu = 0
top += N * nu * sigma
bottom += N * sigma
return ufloat(*(top / bottom))
I'd recommend you read through this recent thread on our discussion forum. The gist is that the definition of keff (and by extension kinf) that we use in OpenMC does not include production of neutrons from (n,xn) in the numerator. So, for your case that would mean that ν=0 when j=(n,xn). "Absorption" in the denominator actually ends up being reduced by (n,2n). More specifically the denominator would be (absorption - (n,2n) - 2*(n,3n) - 3*(n,4n)).
The link you posted above just brings me to a U of I shibboleth. Can you post the DOI instead?
I updated the link to the paper