openmc icon indicating copy to clipboard operation
openmc copied to clipboard

`k_inf` estimator for `IndependentOperator`

Open yardasol opened this issue 2 years ago • 2 comments

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))

yardasol avatar Jul 21 '22 21:07 yardasol

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?

paulromano avatar Jul 21 '22 21:07 paulromano

I updated the link to the paper

yardasol avatar Jul 21 '22 22:07 yardasol