Multigroup discretisation¶

We wish to prove the multigroup formm of the neutron transport equation:

\[\begin{split}\Omega \cdot \nabla \psi_{g}(\mathbf{r}, \Omega) + &\sigma_{g}(\mathbf{r}, \Omega) \psi_{ g}(\mathbf{r}, \Omega) \\ &= \sum_{g'} \sum_{lm} Y_{lm}(\Omega) \sigma^s_{lm, g, g' }(\mathbf{r}) \int_{S^2} Y_{lm}(\Omega') \psi_{g'}(\mathbf{r}, \Omega') d\Omega' + Q_{gp}(\mathbf{r}, \Omega)\end{split}\]

We will show this for the general case (where multiple unknowns per group are allowed) and then specialize it to the standard multigroup transport equation.

Proof¶

We define energy group \(g\) as the energy interval \(E \in [E_{g}, E_{g-1}]\), where \(E_0\) is the highest energy and \(E_G\) is the lowest energy (for a total of \(G\) energy groups).

We now modify the neutron transport equation

\[\begin{split}\Omega \cdot \nabla \psi(\mathbf{r}, \Omega, E) + &\sigma(\mathbf{r}, \Omega) \psi(\mathbf{r}, \Omega, E) \\ &=\iint \sigma^s(\mathbf{r}, \Omega' \to \Omega, E' \to E) \psi(\mathbf{r}, \Omega', E') dE' d\Omega' + Q(\mathbf{r}, \Omega, E)\end{split}\]

as follows. For each group, we define basis functions \(\chi_{gp}\). These basis functions have \(\sum_{p} \chi_{gp} = 1\) for all \(E \in [E_{g}, E_{g-1}]\). Multiply the transport equation by \(\chi_{gp}(E)\) and integrate over \(E\).

Left hand side¶

For the left-hand-side, we have

\[\begin{split}\Omega \cdot \nabla \int_{g} \chi_{gp} \psi(\mathbf{r}, \Omega, E) dE + &\int_{g} \chi_{gp} \sigma(\mathbf{r}, \Omega, E) \psi(\mathbf{r}, \Omega, E) dE \\\end{split}\]

By defining the quantities

\[\begin{split}&\psi_{gp}(\mathbf{r}, \Omega) = \int_{g} \chi_{gp} \psi dE\\ &\sigma_{gp}(\mathbf{r}, \Omega) = \frac{\int_{g} \chi_{gp} \sigma \psi dE}{\psi_{gp}(\mathbf{r}, \Omega)}\end{split}\]

this can be rewritten as

\[\Omega \cdot \nabla \psi_{gp}(\mathbf{r}, \Omega) + \sigma_{gp}(\mathbf{r}, \Omega; \psi) \psi_{gp}(\mathbf{r}, \Omega)\]

Scattering term¶

For the right-hand-side scattering term, we recognize that, in most cases, scattering only depends on \(\Omega'\cdot \Omega\) instead of both separately. Then, we can expand it into a Legendre series:

\[\sigma^s(\mathbf{r}, \Omega' \to \Omega, E' \to E) \equiv \sigma^s(\mathbf{r}, \Omega' \cdot \Omega, E' \to E) =\sum_{l=0}^{\infty} \frac{2l+1}{4\pi} P_l(\Omega \cdot \Omega') \sigma^s_{l}(\mathbf{r}, E' \to E)\]

where \(P_l\) are the Legendre polynomials and \(\sigma^s_{l}(\mathbf{r}, E' \to E)\) are the Legendre moments of the scattering cross section, computed by:

\[\sigma^s_{l}(\mathbf{r}, E' \to E) = 2\pi \int_{-1}^{1} P_l(\mu) \sigma^s(\mathbf{r}, \mu, E' \to E) d\mu\]

This significantly reduces the amount of data required to represent the scattering cross section. Furthermore, this is in practice a good approximation: truncating the sum to the first three terms is already sufficient in many cases. By the spherical harmonic addition theorem and using the definition of the spherical harmonics, this can be rewritten as

\[\sigma^s(\mathbf{r}, \Omega' \to \Omega, E' \to E) = \sum_{l=0}^{\infty} \sum_{m=-l}^{l} Y_{lm}(\Omega) Y_{lm}(\Omega') \sigma^s_{l}(\mathbf{r}, E' \to E)\]

Then, we can fill this into the scattering term (we separate the \(E\) integral into group integrals):

\[\begin{split}&\int_{g} \sum_{g'}\int_{g'} \int_{S^2} \sum_{l=0}^{\infty} \sum_{m=-l}^{l} Y_{lm}(\Omega) Y_{lm}(\Omega') \sigma^s_{l}(\mathbf{r}, E' \to E)\psi(\mathbf{r}, \Omega', E') d\Omega' dE' \chi_{gp}(E) dE \\ &= \sum_{g'} \sum_{lm} Y_{lm}(\Omega) \int_g \int_{g'}\int_{S^2} \sigma^s_{l}(\mathbf{r}, E' \to E) Y_{lm}(\Omega') \psi(\mathbf{r}, \Omega', E') \chi_{gp}(E) dE d\Omega' dE' \\\end{split}\]

Under the integral, we are free to multiply with \(\sum_{g'p'} \chi_{g'p'}(E')=1\). Then, taking the sum out of the integral:

\[\begin{split}\sum_{g'p'} \sum_{lm} Y_{lm}(\Omega) \int_g \int_{g'}\int_{S^2}\chi_{g'p'} (E')\sigma^s_{l}(\mathbf{r}, E' \to E) Y_{lm}(\Omega') \psi(\mathbf{r}, \Omega', E') \chi_{gp}(E) dE d\Omega' dE' \\\end{split}\]

Then, we define:

\[\sigma^s_{lm, gp, g'p'}(\mathbf{r}; \psi) = \frac{\int_{S^2}\int_g \int_{g'}Y_{lm}(\Omega') \sigma^s_{l}(\mathbf{r}, E' \to E) \chi_{gp}(E) \chi_{g'p'}(E') \psi(\mathbf{r}, \Omega', E') dE dE' d\Omega'}{\int_{S^2}\int_{g'} Y_{lm}(\Omega')\psi(\mathbf{r}, \Omega', E') dE'd\Omega' }\]

and we obtain:

\[\sum_{g'p'} \sum_{lm} Y_{lm}(\Omega) \sigma^s_{lm, gp, g'p'}(\mathbf{r}; \psi) \int_{S^2} Y_{lm}(\Omega') \psi_{g'p'}(\mathbf{r}, \Omega') d\Omega'\]

External Source term¶

Finally, for the external source term we can simply define:

\[Q_{gp}(\mathbf{r}, \Omega) = \int_{g} \chi_{gp} Q(\mathbf{r}, \Omega, E) dE\]

Putting it all together¶

This yields the multigroup transport equation, which is a nonlinear equation due to the dependence of the cross sections on the fluxes:

\[\begin{split}\Omega \cdot \nabla \psi_{gp}(\mathbf{r}, \Omega) + &\sigma_{gp}(\mathbf{r}, \Omega; \psi) \psi_{gp}(\mathbf{r}, \Omega) \\ &= \sum_{g'} \sum_{lm} Y_{lm}(\Omega) \sigma^s_{lm , gp, g'p'}(\mathbf{r}; \psi) \int_{S^2} Y_{lm}(\Omega') \psi_{g'p'}(\mathbf{r}, \Omega') d\Omega' + Q_{gp}(\mathbf{r}, \Omega)\end{split}\]

The way out is to use a closure relation. This closure relation approximates the dependence of the cross sections on the fluxes by a good guess. If this guess is exactly right, this consistently yields the exact solution. In practice, one uses cross sections generated from a different code, possibly in a simplified scenario (e.g. infinite medium, homogeneous material compositions, etc.). If the guess is good, then the closure approximation is good, and the result of this equation will be good.

This yields the standard linear multigroup transport equation used in practice:

\[\begin{split}\Omega \cdot \nabla \psi_{gp}(\mathbf{r}, \Omega) + &\sigma_{gp}(\mathbf{r}, \Omega) \psi_{gp}(\mathbf{r}, \Omega) \\ &= \sum_{g'} \sum_{lm} Y_{lm}(\Omega) \sigma^s_{lm, gp, g'p' }(\mathbf{r}) \int_{S^2} Y_{lm}(\Omega') \psi_{g'p'}(\mathbf{r}, \Omega') d\Omega' + Q_{gp}(\mathbf{r}, \Omega)\end{split}\]

Now, there is no dependence of the cross sections on the fluxes anymore. This can be self-consistenly solved by providing appropriate cross sections.

Simplification to single unknown per group¶

If we assume that there is only a single unknown per group (i.e., \(p=1\)), then we can drop the \(p\) index everywhere. This yields the standard multigroup transport equation as used in many textbooks:

with the multigroup cross sections defined as:

\[\begin{split}\sigma_{g}(\mathbf{r}, \Omega) &= \frac{\int_{g} \sigma(\mathbf{r}, \Omega, E) \psi(\mathbf{r}, \Omega, E) dE}{\int_{g} \psi(\mathbf{r}, \Omega, E) dE}\\ \sigma^s_{lm, g, g' }(\mathbf{r}) &= \frac{\int_{S^2}\int_g \int_{g'} Y_{lm}(\Omega')\sigma^s_{l } (\mathbf{r}, E' \to E) \psi(\mathbf{r}, \Omega', E') dE dE' d\Omega'}{\int_{S^2}\int_{g'} Y_{lm}(\Omega')\psi(\mathbf{r}, \Omega', E') dE'd\Omega' }\end{split}\]

It is somewhat surpsising that the total cross sections depends on the angle \(\Omega\). This is because the spectrum is not constant in all directions: then, the flux weighing will yield different results for different angles. In practice, the total cross section is implemented as independent of angle. A correction can be applied (see Lewis & Miller, pg 76) in the scattering operator.

Furthermore, the scattering cross section is now a function of both \(l\) and \(m\) for a similar reason. Physically, this means that neutrons propagating in different directions will have different scattering cross sections. This is unphysical in the first expansion (the material doesn’t care which direction the neutron comes from except for crystals and other anisotropic materials), however, here it arises due to the flux weighing again. In practice, this is not always taken into account and all \(m\) moments are set equal.

Finally, the OpenMC tallies for the scattering cross sections are different: they do not weigh with \(Y_{lm}\) at all. This is somewhat confusing to me, as even in the 1D case in Lewis & Miller eq. 2-57, the cross section is weighed by \(\phi_l\) (and thus with \(P_l\)). My guess is that if the flux profile is sufficiently flat, (\(\phi_{lm}(E)\approx \phi_{lm}\) within the energy group) then both result in the same result, and thus the OpenMC tallies are a good approximation.