1. Background on the Linear Boltzmann Equation

The following textbooks are a good source of information on transport theory Duderstadt and Martin [DM79], computational methods for transport Lewis and Miller [LM84], and nuclear reactor theory Bell and Glasstone [BG79], Duderstadt and Hamilton [DH76]. For a review article on transport approximations, consult Sanchez and McCormick [SM82], Azmy and Sartori [AS10].

1.1. Definitions

In radiation transport, particles are described by their position in space and by their momentum vector. This leads to a six-dimensional phase-space. This increased dimensionality is what makes radiation transport problems computationally expensive.

The six-dimensional phase-space is often described in terms of the particle’s position \(\vec{r}\), the particle’s direction of flight using the unit direction vector \(\vec{\Omega}\), and the particle’s energy \(E\). Hence, to describe the time-dependent distribution of particles in phase-space, we introduce their phase-space density as:

\[n(\vec{r},\vec{\Omega},E,t) \,.\]

It is customary to introduce the following quantities:

• The angular flux \(\psi\):

  \[\psi(\vec{r},\vec{\Omega},E,t) = v(E) \, n(\vec{r},\vec{\Omega},E,t) \,,\]

  where \(v(E)\) is the magnitude of the particle’s velocity.

• The angular current \(\vec{j}\):

  \[\vec{j}(\vec{r},\vec{\Omega},E,t) = \vec{\Omega} \, \psi(\vec{r},\vec{\Omega},E,t) \,.\]

The unit direction vector in \(xyz\) in Cartesian coordinate system is

\[\begin{split}\vec{\Omega} = \begin{bmatrix} \sin{\theta} \cos{\varphi} \\ \sin{\theta} \sin{\varphi} \\ \cos{\theta} \\ \end{bmatrix}\end{split}\]

where \(\theta\) is the polar angle and \(\varphi\) the azimuthal angle. We often use \(\mu = \cos{\theta}\).

1.2. The Linear Boltzmann Equation

The linear Boltzmann transport equation for neutral particles is a conservation statement over the phase-space, stating that the rate of change of the particle’s density is equal to gain terms minus loss terms. It is typically written as follows:

\[\begin{split}\begin{gathered} \frac{1}{v(E)} \frac{\partial \psi(\vec{r},\vec{\Omega},E,t)}{\partial t} = \\ \int dE' \int_{4\pi}d\Omega' \, \sigma_s(\vec{r},\vec{\Omega}'\to \vec{\Omega},E'\to E,t)\psi(\vec{r},\vec{\Omega}',E',t) + Q_{\text{ext}}(\vec{r},\vec{\Omega},E,t) \\ - \vec{\Omega} \cdot \vec{\nabla} \psi(\vec{r},\vec{\Omega},E,t) - \sigma_t(\vec{r},E,t)\psi(\vec{r},\vec{\Omega},E,t) \end{gathered}\end{split}\]

for \(\vec{r} \in \mathcal{D}\) (the spatial domain), \(E\in\mathcal{E}=[0,E_{\text{max}}]\) (the energy range), and \(\vec{\Omega}\in \mathcal{S}^2\) (the unit sphere).

• \(Q_{\text{ext}}\) represents the external source rate density,

• \(\sigma_t(\vec{r},E,t)\) is the total interaction cross section,

• \(\sigma_s(\vec{r},\vec{\Omega}'\to \vec{\Omega},E'\to E,t)\) is the double differential scattering cross section; note that additional energy-angle redistribution events can be modeled, such as \((\text{n},\text{xn})\) neutron interactions; finally, note that for isotropic materials (which is an assumption we make), the angular distribution of the outgoing particles only depends on the cosine of the incoming and outgoing directions, \(\mu_0=\vec{\Omega}'\cdot \vec{\Omega}\), so we will replace \(\sigma_s(\vec{r},\vec{\Omega}'\to \vec{\Omega},E'\to E,t)\) with \(\sigma_s(\vec{r},\vec{\Omega}'\cdot \vec{\Omega},E'\to E,t)\).

In succinct operator notation, the above equation can be re-written as

\[\frac{1}{v} \frac{\partial \Psi}{\partial t} = H\Psi + Q_{\text{ext}} - L\Psi \,,\]

where \(L\) is the streaming + interaction operator and \(H\) is the collision operator.

If fission interactions and delayed neutron production are to be included, the above equation is amended as follows

\[\frac{1}{v} \frac{\partial \Psi}{\partial t} = H\Psi +P_p\Psi + Q_{\text{ext}} + \sum_i S_{d,i} - L\Psi \,,\]

where the prompt fission operator is

\[P_p\Psi = \frac{\chi_p(E)}{4\pi} \int dE \int_{4\pi}d\Omega' \, \nu_p\sigma_f(\vec{r},E'\to E,t)\psi(\vec{r},\vec{\Omega}',E',t) \,,\]

and the delayed neutron operator is

\[S_{d,i} = \frac{\chi_{d,i}(E)}{4\pi} C_i(\vec{r},t) \,,\]

where \(C_i\), the neutron precursor concentration in delayed group \(i\), satisfies its own conservation law.

1.3. Boundary Conditions

Boundary conditions are usually supplied either as a known incoming flux (including a zero function for vacuum boundaries) or as an albedo condition. We denote by \(\Gamma\) the boundary of the spatial domain \(\mathcal{D}\) (that is, \(\Gamma = \partial \mathcal{D}\)).

The following types of boundary conditions are employed:

• The imposed-flux conditions can be expressed as

  \[\psi(\vec{r},\vec{\Omega},E,t) = \psi_{\text{inc}}(\vec{r},\vec{\Omega},E,t) \qquad \forall \vec{r} \in \Gamma^-\]

  where

  \[\Gamma^- = \big\{ \vec{r} \in \Gamma \text{ such that } \vec{\Omega}\cdot\vec{n}(\vec{r}) < 0 \big\}\]

  with \(\vec{n}(\vec{r})\) the outward pointing unit normal vector at position \(\vec{r}\).

• The albedo condition can be expressed as

  \[\psi(\vec{r},\vec{\Omega},E,t) = \mathcal{A}(\vec{r},\vec{\Omega}'\to\vec{\Omega},E'\to E,t) \psi(\vec{r},\vec{\Omega}',E',t)\]

  with

  • \(\mathcal{A}(\vec{r},\vec{\Omega}'\to\vec{\Omega},E'\to E,t)\) the albedo operator, and

  • \(\vec{\Omega}'\cdot \vec{n}(\vec{r}) >0\) (outgoing direction) and the incoming direction is \(\vec{\Omega}=\vec{\Omega}'-2\left(\vec{\Omega}\cdot\vec{n}(\vec{r})\right) \vec{n}(\vec{r})\).

For \(k\)-eigenvalue problems, the boundary conditions usually devolve to a zero-incoming flux (\(\psi_{\text{inc}}=0\)) or a unity albedo (i.e., symmetry line, with \(\mathcal{A}=1\)).

1.4. Initial Conditions

For time-dependent problems, initial conditions are supplied as

\[\psi(\vec{r},\vec{\Omega},E,t=0) = f_0(\vec{r},\vec{\Omega},E) \qquad \forall \vec{r}\in \mathcal{D},\ \forall E \in \mathcal{E}, \ \forall\vec{\Omega}\in \mathcal{S}^2\]

1.5. Expansion of the Angle Redistribution Term

The angle redistribution term

\[\int_{4\pi}d\Omega' \, \sigma_s(\vec{r},\vec{\Omega}'\cdot \vec{\Omega},E'\to E,t)\psi(\vec{r},\vec{\Omega}',E',t)\]

can be expanded in angle on the unit sphere \(\mathcal{S}^2\) using (real-valued) spherical harmonic functions to yield

\[\sum_{\ell=0}^{L_{\text{max}}} \sum_{m=-\ell}^{m=\ell} \, \frac{2\ell+1}{4\pi}\sigma_{s,\ell}(\vec{r},E'\to E,t) Y_{\ell,m}(\vec{\Omega}) \phi_{\ell,m}(\vec{r},E',t)\]

where

• \(L_{\text{max}}\) is the highest order of scattering anisotropy retrained in the cross section expansion (usually, a user-supplied value)

• the Legendre moments of the scattering cross section are defined as

  \[\sigma_{s,\ell}(\vec{r},E'\to E,t) = 2\pi \int_{-1}^1 d\mu_0 \, \sigma_s(\vec{r},\mu_0,E'\to E,t) P_\ell(\mu_0)\]

  where \(\mu_0=\vec{\Omega}'\cdot \vec{\Omega}\) and \(P_\ell(\mu)\) is the Legendre polynomial of degree \(\ell\).

• the moment of the angular flux are defined as

  \[\phi_{\ell,m}(\vec{r},E,t) = \int_{4\pi} d\Omega \, Y_{\ell,m}(\vec{\Omega})\psi(\vec{r},\vec{\Omega},E,t)\]

• the real-valued spherical harmonics are

  \[\begin{split}Y_{\ell,m}(\vec{\Omega}) = \begin{cases} (-1)^m \sqrt(2)\sqrt{ \frac{(2\ell + 1)}{4\pi} \frac{(\ell-|m|)!}{(\ell+|m|)!}}P_{\ell}^{|m|}(\cos\theta)\sin{|m|\varphi} & \text{if } m < 0 \\ \\ \sqrt{ \frac{(2\ell + 1)}{4\pi}} P_{\ell}^{m}(\cos\theta) & \text{if } m = 0 \\ \\ (-1)^m \sqrt(2)\sqrt{ \frac{(2\ell + 1)}{4\pi} \frac{(\ell-m)!}{(\ell+m)!}}P_{\ell}^{m}(\cos\theta)\cos{m\varphi} & \text{if } m > 0 \\ \end{cases}\end{split}\]

  where \(P^m_\ell(\mu)\) is the associated Legendre function of degree \(\ell\) and of order \(m\).

In operator notation, the steady-state, source-driven problem \(L\Psi = H\Psi + Q_{\text{ext}}\) can now be re-written as

\[L\Psi = M\Sigma \Phi + Q_{\text{ext}} \qquad \text{with } \Phi = D \Psi\]

where

• \(\Phi\) are the flux moments,

• \(D\) is the discrete-to-moment operator and denotes the angular integration of the angular flux to yield the flux moment (\(\Phi = D \Psi\)). The term discrete stems from the fact that the integration is performed using a quadrature rule, hence at discrete directions for the angular flux,

• \(\Sigma\) denotes the matrix containing the Legendre moments of the scattering cross section, and

• \(M\) denotes the moment-to-discrete operator, which takes the source moments (\(\Sigma \Phi\)) and evaluates that term in direction.

If we denote by \(N_{\text{mom}}\) the maximum number of flux moments, we have

\[\begin{split}\begin{aligned} N_{\text{mom}} &= L_{\text{max}}+1 &\text{in 1D,} \\ N_{\text{mom}} &= \frac{(L_{\text{max}}+1)(L_{\text{max}}+2)}{2} &\text{in 2D,} \\ N_{\text{mom}} &= (L_{\text{max}}+1)^2 &\text{in 3D.} \end{aligned}\end{split}\]

With the introduction of the moment-to-discrete operator, we can also update the fission production operator

\[\text{from}\quad P\Psi \quad\text{to}\quad M F \Phi\]

where, hereafter, \(F\) will denote the fission operator acting on flux moments.

1.6. Short summary of Transport Equations Solved in OpenSn

• Time-dependent transport problem:

  \[\frac{1}{v} \frac{\partial \Psi}{\partial t} = M \Sigma \Phi + MF_p\Phi + Q_{\text{ext}} + \sum_i S_{d,i} - L\Psi \,,\]

• Steady-state, subcritical, source-driven transport problem:

  \[L\Psi = M \Sigma \Phi + M F \Phi + Q_{\text{ext}} \,,\]

  where \(P\) is the total (prompt+delayed) fission production operator.

• \(k\)-eigenvalue transport problem:

  \[L\Psi = M \Sigma \Phi + \frac{1}{k_{\text{eff}}} M F \Phi\]

1.7. Streaming Term in Cartesian and Curvilinear Coordinate Systems

In \(xyz\) coordinates, the streaming term is given by

\[\vec{\Omega}\cdot \vec{\nabla}\psi = \Omega_x \partial_x\psi + \Omega_y \partial_y\psi + \Omega_z \partial_z\psi\]

In \(r\theta z\) cylindrical coordinates, the streaming term is given by

\[\vec{\Omega}\cdot \vec{\nabla} = \frac{\xi}{r} \partial_r(r\psi) - \frac{1}{r}\partial_\varphi (\eta\psi) + \mu \partial_z\psi\]

with \(\xi=\vec{\Omega} \cdot \vec{e}_r=\sqrt{1-\mu^2}\cos{\varphi}\), \(\eta=\vec{\Omega} \cdot \vec{e}_\theta=\sqrt{1-\mu^2}\sin{\varphi}\), \(\mu=\vec{\Omega}\cdot \vec{e}_z\)

1.8. References

[DM79]

James J Duderstadt and William Russell Martin. Transport theory. John Wiley & Sons, 1979.

[LM84]

Elmer Eugene Lewis and Warren F Miller. Computational methods of neutron transport. John Wiley and Sons, Inc., New York, NY, 1984. ISBN 0-89448-452-4.

[BG79]

George I Bell and Samuel Glasstone. Nuclear reactor theory. RE Krieger Publishing Company, 1979. ISBN 978-0-442-20684-0.

[DH76]

James J Duderstadt and Louis J Hamilton. Nuclear reactor analysis. Wiley, 1976.

[SM82]

R Sanchez and Norman J McCormick. Review of neutron transport approximations. Nuclear Science and Engineering, 1982.

[AS10]

Y. Azmy and E. Sartori. Nuclear Computational Science: A Century in Review. Springer Netherlands, 2010. ISBN 9789048134113. URL: https://books.google.com/books?id=4imnyVPQ-0AC.