2. Method of Characteristics¶

The method of characteristics (MOC) is a widely used technique for solving partial differential equations, including the Boltzmann form of the neutron transport equation [Askew]. MOC is used to solve the transport equation in 2D by discretizing both polar and azimuthal angles and integrating the multi-group characteristic form of the equation for a particular azimuthal and polar angle quadrature. The following sections detail the derivation of the characteristic form of the 2D neutron transport equation solved in the OpenMOC method of characteristics scheme [Boyd].

Section 2.1 introduces the Boltzmann form of the neutron transport equation parametrized in 6-dimensional phase space over position, angle and energy. The following several sections introduce the various approximations made to this equation:

characteristic transformation (Section 2.2)
energy discretization (Section 2.3)
discrete ordinates approximation (Section 2.4)
isotropic scattering approximation (Section 2.5)
flat source approximation (Section 2.6)
constant cross-section approximation (Section 2.7)
integrating factor solution (Section 2.8)
track area approximation (Section 2.9)
azimuthal planar projection (Section 2.10)

The final equations applied in OpenMOC to solve for the FSR source and scalar flux derived in the following sections are summarized below:

The source in each flat source region

$\boxed{Q_{i,g} = \frac{1}{4\pi}\left(\displaystyle\sum\limits_{g'=1}^G \Sigma^S_{i,g'\rightarrow g}\Phi_{i,g'} + \frac{\chi_{i,g}}{k_{eff}}\displaystyle\sum\limits_{g'=1}^G\nu\Sigma^F_{i,g'}\Phi_{i,g'}\right)}$

Change in angular flux along a track segment

$\boxed{\Delta\Psi_{k,i,g,p} = \Psi_{k,g,p}(s') - \Psi_{k,g,p}(s'') = \left(\Psi_{k,g,p}(s') - \frac{Q_{i,g}}{\Sigma^T_{i,g}}\right)(1 - e^{-\tau_{k,i,g,p}})}$

Scalar flux in each flat source region

$\boxed{\Phi_{i,g} = \frac{4\pi}{\Sigma_{i,g}}\left[Q_{i,g} + \frac{1}{A_i}\displaystyle\sum\limits_{k\in A_{i}}\displaystyle\sum\limits_{p=1}^{P}\omega_{m(k)}\omega_{p}\omega_{k}\sin\theta_{p}\Delta\Psi_{k,i,g,p}\right]}$

2.1. Introduction to the Boltzmann Equation¶

The Boltzmann form of the steady-state neutron transport equation is given by the following:

(1) $\mathbf{\Omega} \cdot \nabla \Psi(\mathbf{r},\mathbf{\Omega},E) + \Sigma^T(\mathbf{r},E)\Psi(\mathbf{r},\mathbf{\Omega},E) = \int_{0}^{\infty} \mathrm{d}E' \int_{4\pi} \mathrm{d}\mathbf{\Omega'}\Sigma^S(\mathbf{r},{\mathbf{\Omega'}\rightarrow\mathbf{\Omega}},{E'\rightarrow E}) \Psi(\mathbf{r},\mathbf{\Omega'},E') + \frac{\chi(\mathbf{r},E)}{4\pi k_{eff}} \int_{0}^{\infty} \mathrm{d}E' \nu\Sigma^F(\mathbf{r},E') \int_{4\pi} \mathrm{d}\mathbf{\Omega'}\Psi(\mathbf{r},\mathbf{\Omega'},E')$

Each of the variables in use is defined in Table 1. This is a balance equation between neutrons lost to transport, lost to absorption, produced or lost from scattering and those produced from fission. It should be noted that this equation assumes isotropic emission from fission.

Variable	Description
$\mathbf{r}$	Spatial position vector
$\mathbf{\Omega}$	Angular direction vector
$E$	Neutron energy
$\Psi$	Angular neutron flux
$k_{eff}$	Effective neutron multiplication factor
$\Sigma^T$	Neutron total cross-section
$\Sigma^S$	Neutron scattering cross-section
$\Sigma^F$	Neutron fission cross-section
$\chi$	Energy spectrum for fission neutrons
$\nu$	Average number of neutrons emitted per fission

Table 1: Variables in the Boltzmann equation.

The first step is to simplify this equation by defining those quantities on the right hand side as the total neutron source $Q(\mathbf{r},\mathbf{\Omega},E)$ :

(2) $Q(\mathbf{r},\mathbf{\Omega},E) = \int_{0}^{\infty} \mathrm{d}E' \int_{4\pi} \mathrm{d}\mathbf{\Omega'}\Sigma^S(\mathbf{r},{\mathbf{\Omega'}\rightarrow\mathbf{\Omega}},{E'\rightarrow E}) \Psi(\mathbf{r},\mathbf{\Omega'},E') + \frac{\chi(\mathbf{r},E)}{4\pi k_{eff}} \int_{0}^{\infty} \mathrm{d}E' \int_{4\pi} \mathrm{d}\mathbf{\Omega'} \nu\Sigma^F(\mathbf{r},E')\Psi(\mathbf{r},\mathbf{\Omega'},E')$

The transport equation can now be more concisely written as follows:

(3) $\mathbf{\Omega} \cdot \nabla \Psi(\mathbf{r},\mathbf{\Omega},E) + \Sigma^T(\mathbf{r},E)\Psi(\mathbf{r},\mathbf{\Omega},E) = Q(\mathbf{r},\mathbf{\Omega},E)$

2.2. The Characteristic Transformation¶

The characteristic form of the Boltzmann equation is found by a change of variables by parametrizing $\mathbf{r}$ with respect to some reference location $\mathbf{r_0}$ :

(4) $\mathbf{r} = (x(s), y(s)) = (x_0+s\mathbf{\Omega_x}, y_0+s\mathbf{\Omega_y}) = \mathbf{r_0}+s\mathbf{\Omega}$

For any location $\mathbf{r}$ of interest, each angular direction vector $\mathbf{\Omega'}$ is matched to a corresponding reference location $\mathbf{r_{0}'}$ defined such that $\mathbf{r} = \mathbf{r_{0}'} + s\mathbf{\Omega'}$ . This parametrization for position may be substituted into the source and transport equations to obtain the following form for each:

(5) $Q(\mathbf{r},\mathbf{\Omega},E) = \int_{0}^{\infty} \mathrm{d}E' \int_{4\pi} \mathrm{d}\mathbf{\Omega'}\Sigma^S(\mathbf{r_0'}+s\mathbf{\Omega'},{\mathbf{\Omega'}\rightarrow\mathbf{\Omega}},{E'\rightarrow E}) \Psi(\mathbf{r_0'}+s\mathbf{\Omega'},\mathbf{\Omega'},E') + \frac{\chi(\mathbf{r_0}+s\mathbf{\Omega},E)}{4\pi k_{eff}} \int_{0}^{\infty} \mathrm{d}E' \nu\Sigma^F(\mathbf{r_0}+s\mathbf{\Omega},E') \int_{4\pi} \mathrm{d}\mathbf{\Omega'} \Psi(\mathbf{r_0'}+s\mathbf{\Omega'},\mathbf{\Omega'},E')$

(6) $\mathbf{\Omega} \cdot \nabla \Psi(\mathbf{r_0}+s\mathbf{\Omega},\mathbf{\Omega},E) + \Sigma^T(\mathbf{r_0}+s\mathbf{\Omega},E)\Psi(\mathbf{r_0}+s\mathbf{\Omega},\mathbf{\Omega},E) = Q(\mathbf{r_0}+s\mathbf{\Omega},\mathbf{\Omega},E)$

Applying the differential operator to the angular flux in (6) leads to the characteristic form of the Boltzmann equation:

(7) $\frac{d}{ds}\Psi(\mathbf{r_0}+s\mathbf{\Omega},\mathbf{\Omega},E) + \Sigma^T(\mathbf{r_0}+s\mathbf{\Omega},E)\Psi(\mathbf{r_0}+s\mathbf{\Omega},\mathbf{\Omega},E) = Q(\mathbf{r_0}+s\mathbf{\Omega},\mathbf{\Omega},E)$

For brevity, the remainder of this section will assume the dependence of $s$ on the reference position $\mathbf{r_0}$ and $\mathbf{\Omega}$ and will simplify this as $\mathbf{r_0} + s\mathbf{\Omega} \rightarrow s$ such that the characteristic equation can be written as the following:

(8) $\frac{d}{ds}\Psi(s,\mathbf{\Omega},E) + \Sigma^T(s,E)\Psi(s,\mathbf{\Omega},E) = Q(s,\mathbf{\Omega},E)$

This equation can be solved through the use of an integrating factor:

(9) $e^{-\int_{0}^s\mathrm{d}s'\Sigma^T(s',E)}$

The final analytical solution to the characteristic equation is therefore:

(10) $\Psi(s,\mathbf{\Omega},E) = \Psi(\mathbf{r_{0}},\mathbf{\Omega},E)e^{-\int_{0}^s\mathrm{d}s'\Sigma^T(s',E)} + \int_0^s\mathrm{d}s''Q(s'',\mathbf{\Omega},E)e^{-\int_{s''}^s\mathrm{d}s'\Sigma^T(s',E)}$

2.3. The Multi-Group Energy Approximation¶

Equation (10) is defined with $\Psi$ , $Q$ and $\Sigma^T$ as continuous functions of energy. The first approximation to numerically solve this equation is to discretize the energy domain into distinct energy groups $g \in G = \{1, 2, ..., G\}$ where group $g$ spans the continuous range of energies from $E_{g}$ to $E_{g-1}$ . This is otherwise known as the multi-group approximation. The multi-group form of the Boltzmann equation is presented below:

(11) $\mathbf{\Omega} \cdot \nabla \Psi_g(s,\mathbf{\mathbf{\Omega}}) + \Sigma^T_{g}(s)\Psi_g(s,\mathbf{\Omega}) = Q_g(s,\mathbf{\Omega})$

The characteristic form of the equation given in (8) can also be written in multi-group form:

(12) $\frac{d}{ds}\Psi_{g}(s,\mathbf{\Omega}) + \Sigma^T_{g}(s)\Psi_{g}(s,\mathbf{\Omega}) = Q_g(s,\mathbf{\Omega})$

Likewise, the multi-group form of the neutron source (5) is given by:

(13) $Q_g(s,\mathbf{\Omega}) = \displaystyle\sum\limits_{g'=1}^G \int_{4\pi} \mathrm{d}\mathbf{\Omega'}\Sigma_{g'\rightarrow g}^S(s,{\mathbf{\Omega'}\rightarrow\mathbf{\Omega}}) \Psi_{g'}(s,\mathbf{\Omega'}) + \frac{\chi_{g}(s)}{4\pi k_{eff}} \displaystyle\sum\limits_{g'=1}^G \nu\Sigma_{g'}^F(s) \int_{4\pi} \mathrm{d}\mathbf{\Omega'} \Psi_{g'}(s,\mathbf{\Omega'})$

It directly follows from (10) and (8) that the solution to the multi-group characteristic neutron transport equation is the following:

(14) $\Psi_g(s,\mathbf{\Omega}) = \Psi_g(\mathbf{r_{0}},\mathbf{\Omega})e^{-\int_{0}^s\mathrm{d}s'\Sigma_g^T(s')} + \int_0^s\mathrm{d}s''Q_g(s'',\mathbf{\Omega})e^{-\int_{s''}^s\mathrm{d}s'\Sigma_g^T(s')}$

Where both (14) and (13) make use of the energy condensed cross-sections $\Sigma^T$ , $\Sigma^F$ , $\Sigma^S$ , and $\chi$ :

(15) $\Sigma_{g}^T(s) = \frac{\int_{E_{g}}^{E_{g-1}}\mathrm{d}E'\Sigma^T(s,E')\Psi(s,\mathbf{\Omega},E')}{\int_{E_{g}}^{E_{g-1}}\mathrm{d}E'\Psi(s,\mathbf{\Omega},E')}$

(16) $\Sigma_{g}^F(s) = \frac{\int_{E_{g}}^{E_{g-1}}\mathrm{d}E'\Sigma^F(s,E')\Psi(s,\mathbf{\Omega},E')}{\int_{E_{g}}^{E_{g-}1}\mathrm{d}E'\Psi(s,\mathbf{\Omega},E')}$

(17) $\Sigma_{g'\rightarrow g}^S(s,\mathbf{\Omega'}\rightarrow \mathbf{\Omega}) = \frac{\int_{E_{g'}}^{E_{g'-1}}\mathrm{d}E'\int_{E_{g}}^{E_{g-1}}\mathrm{d}E''\Sigma^S(s,\mathbf{\Omega'}\rightarrow \mathbf{\Omega},E'\rightarrow E'')\Psi(s,\mathbf{\Omega'},E')}{\int_{E_{g'}}^{E_{g'-1}}\mathrm{d}E'\Psi(s,\mathbf{\Omega'},E')}$

(18) $\chi_{g'\rightarrow g}(s) = \frac{\int_{E_{g'}}^{E_{g'-1}}\mathrm{d}E'\int_{E_{g}}^{E_{g-1}}\mathrm{d}E''\chi(s,E'\rightarrow E'')\nu\Sigma^F(s,\mathbf{\Omega},E')\Psi(s,\mathbf{\Omega'},E')}{\int_{E_{g'}}^{E_{g'-1}}\mathrm{d}E'\nu\Sigma^F(s,\mathbf{\Omega},E')\Psi(s,\mathbf{\Omega'},E')}$

Although (18) assumes a dependence of $\chi$ on both the energy of the neutron causing fission $g'$ and the fission emission energy group $g$ , the former is typically summed over to simplify the multi-group $\chi$ to the following approximation:

(19) $\chi_{g}(s) = \displaystyle\sum\limits_{g=1}^{G}\chi_{g'\rightarrow g}(s)$

2.4. The Discrete Ordinates Approximation¶

The discrete ordinates approximation is introduced to approximate the integral over the angular domain in the source (13). This is equivalent to applying quadrature rules to evaluate the integral over the angular flux using a weighted sum of fluxes at specific angles where weights $w_{m}$ are introduced for each of the quadrature points $\mathbf{\Omega_{m}} \; \forall \; m \in \{1, ..., M\}$ .

(20) $\Phi_{g}(s) = \int_{4\pi}\mathrm{d}\mathbf{\Omega'}\Psi_{g}(s,\mathbf{\Omega'}) \approx \displaystyle\sum\limits_{m=1}^{M}w_{m}\Psi_{g}(s,\mathbf{\Omega_{m}})$

The integrated angular flux $\Phi_{g}(s)$ is termed the scalar flux. Substituting this approximation to the angular flux integral into (13) leads to the following approximation to the source $Q_{m,g}(s) \approx Q_{g}(s,\mathbf{\Omega_{m}})$ at each quadrature point $\mathbf{\Omega_{m}}$ :

(21) $Q_{m,g}(s) = \displaystyle\sum\limits_{g'=1}^G \displaystyle\sum\limits_{m'=1}^{M}w_{m'}\Sigma_{g'\rightarrow g}^S(s,{\mathbf{\Omega_{m'}}\rightarrow\mathbf{\Omega_{m}}}) \Psi_{g'}(s,\mathbf{\Omega_{m'}}) + \frac{\chi_{g}(s)}{4\pi k_{eff}} \displaystyle\sum\limits_{g'=1}^G \displaystyle\sum\limits_{m'=1}^{M}w_{m'}\nu\Sigma_{g'}^F(s)\Psi_{g'}(s,\mathbf{\Omega_{m'}})$

Substituting this approximation to the source into (14) one obtains the characteristic solution for the angular flux $\Psi_{m,g}(s) \approx \Psi_g(s,\mathbf{\Omega_{m}})$ at each quadrature point $\mathbf{\Omega_{m}}$ :

(22) $\Psi_{m,g}(s) = \Psi_{m,g}(\mathbf{r_{0}})e^{-\int_{0}^{s}\mathrm{d}s'\Sigma_g^T(s')} + \int_0^{s_{m}}\mathrm{d}s''Q_{m,g}(s'')e^{-\int_{s''}^{s}\mathrm{d}s'\Sigma_g^T(s')}$

Equations (21) and (22) may be further decomposed into azimuthal and polar angle quadratures $m \in \{1, 2, ..., M\}$ and $p \in \{1, 2, ..., P\}$ with weights $w_{m}$ and $w_{p}$ for the azimuthal plane and axial dimension, respectively:

(23) $Q_{m,p,g}(s) = \displaystyle\sum\limits_{g'=1}^G \displaystyle\sum\limits_{m'=1}^{M} \displaystyle\sum\limits_{p'=1}^{P} w_{m'}w_{p'}\Sigma_{g'\rightarrow g}^S(s,{\mathbf{\Omega_{m',p'}}\rightarrow\mathbf{\Omega_{m,p}}}) \Psi_{g'}(s,\mathbf{\Omega_{m',p'}}) + \frac{\chi_{g}(s)}{4\pi k_{eff}} \displaystyle\sum\limits_{g'=1}^G \displaystyle\sum\limits_{m'=1}^{M} \displaystyle\sum\limits_{p'=1}^{P} w_{m'}w_{p'} \nu\Sigma_{g'}^F(s)\Psi_{g'}(s,\mathbf{\Omega_{m',p'}})$

(24) $\Psi_{m,p,g}(s) = \Psi_{m,p,g}(\mathbf{r_{0}})e^{-\int_{0}^s\mathrm{d}s'\Sigma_g^T(s')} + \int_0^s\mathrm{d}s''Q_{m,p,g}(s'')e^{-\int_{s''}^s\mathrm{d}s'\Sigma_g^T(s')}$

2.5. The Isotropic Scattering Approximation¶

An additional approximation that is made to simplify the evaluation of the source in (23) is to assume that the scattering source is isotropic. This approximation allows the total source to be expressed solely in terms of the scalar flux:

(25) $Q_{g}(s) = \frac{1}{4\pi}\left(\displaystyle\sum\limits_{g'=1}^G \Sigma^S_{g'\rightarrow g}(s)\Phi_{g'}(s) + \frac{\chi_{g}(s)}{k_{eff}}\displaystyle\sum\limits_{g'=1}^G\nu\Sigma^F_{g'}(s)\Phi_{g'}(s)\right)$

The subscripts $m$ and $p$ for the azimuthal and polar angles, respectively, have been dropped from $Q_{g}(s)$ since they have been embedded in the integral over angular phase space to obtain the scalar flux $\Phi_{g}(s)$ .

2.6. The Flat Source Region Approximation¶

Another common approximation for MOC is to assume that the source $Q_g$ is constant across discrete spatial cells termed flat source regions (FSRs). This implies that the source does not vary along a characteristic $k$ entering FSR $i$ at $s'$ and exiting at $s''$ :

(26) $Q_{i,g} = Q_{g}(s') = Q_{g}(s'') = Q_{g}(s) \;\;\; , \;\;\; s \in [s', s'']$

2.7. The Linear Source Region Approximation¶

A more accurate description of the source in spatial regions is to assume a linear variation. This is typically sufficient for the moderator in a PWR, when each channel is also cut in azimuthal source regions. The source then varies along each characteristic lines. The reader should refer themselves to [Ferrer] and [Gunow] for more details on the track-based linear source approximation and its implementation in OpenMOC.

(27) $Q_{g}(s) = q_{t,g,0} + q_{t,g,1} (s - l_{t} / 2)$

$l_{t}$ is the length of the segment considered, while $q_{g,0}$ and $q_{g,1}$ are track dependent coefficients that describe the source.

2.8. The Constant Cross-Section Approximation¶

In addition to the flat source approximation, it is assumed that the material properties are constant across each FSR. The area-averaged cross-sections for FSR $i \in \{1, 2, ..., I\}$ with area $A_{i}$ are defined as:

(28) $\Sigma_{i,g}^{T} = \frac{\int_{\mathbf{r}\in A_{i}}\mathrm{d}\mathbf{r}\Sigma_{g}^T(\mathbf{r})\Phi_{g}(\mathbf{r})}{\int_{\mathbf{r}\in A_{i}}\mathrm{d}\mathbf{r}\Phi_{g}(\mathbf{r})}$

(29) $\Sigma_{i,g}^{F} = \frac{\int_{\mathbf{r}\in A_{i}}\mathrm{d}\mathbf{r}\Sigma_{g}^F(\mathbf{r})\Phi_{g}(\mathbf{r})}{\int_{\mathbf{r}\in A_{i}}\mathrm{d}\mathbf{r}\Phi_{g}(\mathbf{r})}$

(30) $\Sigma_{i,g'\rightarrow g}^{S} = \frac{\int_{\mathbf{r}\in A_{i}}\mathrm{d}\mathbf{r}\Sigma_{g'\rightarrow g}^S(\mathbf{r})\Phi_{g'}(\mathbf{r})}{\int_{\mathbf{r}\in A_{i}}\mathrm{d}\mathbf{r}\Phi_{g'}(\mathbf{r})}$

(31) $\chi_{i,g} = \frac{\int_{\mathbf{r}\in A_{i}}\mathrm{d}\mathbf{r}\chi_{g}(\mathbf{r})}{\int_{\mathbf{r}\in A_{i}}\mathrm{d}\mathbf{r}}$

The flat source term $Q_{i,g}$ for FSR $i$ with area $A_i$ is defined in terms of both fission and scattering from the area-averaged scalar flux $\Phi_{g,i}$ within the FSR:

(32) $Q_{i,g} = \frac{1}{4\pi}\left(\displaystyle\sum\limits_{g'=1}^G \Sigma^S_{i,g'\rightarrow g}\Phi_{i,g'} + \frac{\chi_{i,g}}{k_{eff}}\displaystyle\sum\limits_{g'=1}^G\nu\Sigma^F_{i,g'}\Phi_{i,g'}\right)$

(33) $\Phi_{i,g} = \frac{\int_{\mathbf{r}\in A_{i}}\mathrm{d}\mathbf{r}\Phi_{g}(\mathbf{r})}{\int_{\mathbf{r}\in A_{i}}\mathrm{d}\mathbf{r}}$

The multi-group nuclear cross-sections for each FSR are an input to OpenMOC. As a result, the area-averaging integrals must be performed by some pre-processing method such as Monte Carlo.

2.9. The Integrating Factor Solution¶

Each chracteristic may be discretized into segments across individual FSRs. This approximation allows (24) to be localized to a segment of characteristic $k$ across FSR $i$ from its entry point at $s'$ to exit point at $s''$ . By defining the integrating factor in terms of the optical length $\tau_{k,i,g} = \Sigma^T_{i,g}(s''-s')$ one may analytically evaluate the integrals in (24) and express the outgoing flux along the characteristic as follows:

(34) $\Psi_{k,g}(s'') = \Psi_{k,g}(s')e^{-\tau_{k,i,g}} + \frac{Q_{i,g}}{\Sigma^T_{i,g}}(1 - e^{-\tau_{k,i,g}})$

With minor algebraic rearrangement, the change in angular flux along the characteristic is given by the following:

(35) $\Delta\Psi_{k,g} = \Psi_{k,g}(s') - \Psi_{k,g}(s'') = \left(\Psi_{k,g}(s') - \frac{Q_{i,g}}{\Sigma^T_{i,g}}\right)(1 - e^{-\tau_{k,i,g}})$

2.10. The Track Area Approximation¶

The key quantity remaining to be determined is the integral over area for the FSR area-averaged scalar flux $\Phi_{g,i}$ in (33). The track area approximation is used to compute this value numerically.

First, define $l_{k,i}=s''-s'$ such that the average angular flux in FSR $i$ along characteristic $k$ is the following integral:

(36) $\overline{\Psi}_{k,i,g} = \frac{1}{l_{k,i}}\int_{s'}^{s''} \Psi_{k,i,g}(s) \mathrm{d}s$

Upon evaluating the integral, the average angular flux along the characteristic can be reduced to the following algebraic expression:

(37) $\overline{\Psi}_{k,i,g} = \frac{1}{l_{k,i}}\left[\frac{\Psi_{k,g}(s')}{\Sigma_{i,g}^T}(1 - e^{-\tau_{k,i,g}}) + \frac{l_{k,i}Q_{i,g}}{\Sigma_{i,g}^T}\left(1 - \frac{(1 - e^{-\tau_{k,i,g}})}{\tau_{k,i,g}}\right)\right]$

Assuming a constant source and cross-sections in FSR $i$ , the value given for the average angular flux in (37) is exact. In order to exactly compute the area-averaged scalar flux, the average angular flux from every characteristic crossing FSR $i$ must be taken into account. This is numerically intractable; hence, an appropriate subset $K$ of characteristics, henceforth known as tracks, is chosen and the integral over the area of the FSR is performed using quadrature rules with a weight $w_{k}$ for each track $k \in K$ crossing through the FSR $k \in A_{i}$ . The contribution $\overline{\Psi}_{k,i,g}$ of track $k$ with azimuthal and polar quadrature weights denoted by $w_{m(k)}$ and $w_{p(k)}$ , respectively, is then integrated to find the area-averaged scalar flux in FSR $i$ as follows:

(38) $\Phi_{i,g} = \frac{\int_{\mathbf{r}\in A_{i}}\mathrm{d}\mathbf{r}\int_{4\pi}\mathrm{d}\mathbf{\Omega}\Psi_{g}(\mathbf{r},\mathbf{\Omega})}{\int_{\mathbf{r}\in A_{i}}\mathrm{d}\mathbf{r}} \approx \frac{4\pi\displaystyle\sum\limits_{k\in A_{i}}w_{m(k)}w_{p(k)}w_{k}l_{k,i}\sin\theta_{p(k)}\overline{\Psi}_{k,i,g}}{\displaystyle\sum\limits_{k\in A_{i}}w_kl_{k,i}\sin\theta_{p(k)}}$

In (38), the angle $\theta_{p(k)}$ formed by characteristic $k$ with respect to the polar axis is introduced to project the length of the characteristic segment $l_{k,i}$ onto the azimuthal plane. In this application of quadrature to approximate an area integral, the weights can be thought of as the effective width of each track $k$ .

The denominator in (38) then simplifies to the area $A_i$ :

(39) $\Phi_{i,g} \approx \frac{4\pi}{A_{i}}\displaystyle\sum\limits_{k\in A_{i}}w_{m(k)}w_{p(k)}w_{k}l_{k,i}\sin\theta_{p(k)}\overline{\Psi}_{k,i,g}$

The scalar flux can be found in terms of average angular fluxes from each track by substituting the expression for the average angular flux from (37) into (39) and rearranging:

(40) $\Phi_{i,g} = \frac{4\pi}{\Sigma_{i,g}}\left[Q_{i,g} + \frac{1}{A_i}\displaystyle\sum\limits_{k\in A_{i}}\omega_{m(k)}\omega_{p(k)}\omega_{k}\sin\theta_{p(k)}\left(\Psi_{k,i,g}(s') - \frac{Q_{i,g}}{\Sigma_{i,g}^T}\right)(1 - e^{-\tau_{k,i,g}})\right]$

The final form for the scalar flux can be simplified in terms of the change in angular flux $\Delta\Psi_{k,i,g}$ along each track segment as defined in (35):

(41) $\Phi_{i,g} = \frac{4\pi}{\Sigma_{i,g}}\left[Q_{i,g} + \frac{1}{A_i}\displaystyle\sum\limits_{k\in A_{i}}\omega_{m(k)}\omega_{p(k)}\omega_{k}\sin\theta_{p(k)}\Delta\Psi_{k,i,g}\right]$

2.11. Projection from the Azimuthal Plane¶

The preceding sections used track segment lengths $l_{k,i}$ in 3D. In practice, the memory footprint for storing track segment data is greatly reduced if the polar angle quadrature is replicated for each azimuthal quadrature point. Such a quadrature allows for track segments to be stored in the 2D azimuthal plane and projected into 3D for each polar angle when necessary. The projection results in some minor changes to the equations presented in the previous sections.

In what follows, each track segment length $l_{k,i}$ will be assumed to reside within the azimuthal plane. Likewise, the optical length $\tau_{k,i,g} = \Sigma^T_{k,i,g}l_{k,i}$ also resides in the azimuthal plane. For notational simplicity, the 3D projection of the track segment length for polar angle $p$ will be denoted by $l_{k,i,p} = \frac{l_{k,i}}{\sin\theta_{p}}$ and the optical length by $\tau_{k,i,g,p} = \Sigma^T_{k,i,g}l_{k,i,p}$ .

First, the polar angle must be accounted for in the expression for the track segment average angular flux to project the segment length into the polar dimension:

(42) $\overline{\Psi}_{k,i,g,p} = \frac{1}{l_{k,i,p}}\left[\frac{\Psi_{k,g,p}(s')}{\Sigma_{i,g}^T}(1 - \exp(-\tau_{k,i,g,p})) + \frac{l_{k,i,p}Q_{i,g}}{\Sigma_{i,g}^T}\left(1 - \frac{(1 - \exp(-\tau_{k,i,g,p}))}{\tau_{k,i,g,p}}\right)\right]$

Next, $\sin\theta_{p(k)}$ is dropped and a summation over polar angles is incorporated into the area-averaged scalar flux in (39):

(43) $\Phi_{i,g} = \frac{4\pi}{A_i}\displaystyle\sum\limits_{k \in A_i}\displaystyle\sum\limits_{p=1}^{P}\omega_{m(k)}\omega_{p}\omega_{k}l_{k,i}\overline{\Psi}_{k,i,g,p}$

The scalar flux can be found in terms of average angular fluxes from each track by substituting the expression for the average angular flux from (42) into (43) and rearranging:

(44) $\Phi_{i,g} = \frac{4\pi}{\Sigma_{i,g}}\left[Q_{i,g} + \frac{1}{A_i}\displaystyle\sum\limits_{k\in A_{i}}\displaystyle\sum\limits_{p=1}^{P}\omega_{m(k)}\omega_{p}\omega_{k}\sin\theta_{p}\left(\Psi_{k,i,g,p}(s') - \frac{Q_{i,g}}{\Sigma_{i,g}^T}\right)(1 - e^{-\tau_{k,i,g,p}})\right]$

The final form for the scalar flux can be simplified in terms of the change in angular flux $\Delta\Psi_{k,i,g,p}$ along each track segment as defined in (41):

(45) $\Phi_{i,g} = \frac{4\pi}{\Sigma_{i,g}}\left[Q_{i,g} + \frac{1}{A_i}\displaystyle\sum\limits_{k\in A_{i}}\displaystyle\sum\limits_{p=1}^{P}\omega_{m(k)}\omega_{p}\omega_{k}\sin\theta_{p}\Delta\Psi_{k,i,g,p}\right]$

This is the form of the transport equation solved by the MOC formulation used in OpenMOC.

2.12. References¶

[Askew]

Askew, “A Characteristics Formulation of the Neutron Transport Equation in Complicated Geometries.” Technical Report AAEW-M 1108, UK Atomic Energy Establishment (1972).

[Boyd]

Boyd, “Massively Parallel Algorithms for Method of Characteristics Neutral Particle Transport on Shared Memory Computer Architectures.” M.S. Thesis, Massachusetts Institute of Technology (2014).

[Ferrer]

Ferrer and J. Rhodes, “A Linear Source Approximation Scheme for the Method of Characteristics,” volume 77, p. 119–136, 1981.

[Gunow]

Gunow “Full Core 3D Neutron Transport Simulation Using the Method of Characteristics with Linear Sources”, PhD Thesis, Massachusetts Institute of Technology (2018).