2.3. One- and Two-Dimensional Spatial Discretization#
A standard sweeping method [BG70] forms the core of the inner iteration discussed in the multigroup discrete ordinates method. For this process, the angular flux is calculated at the cell boundary and used to average the flux at the cell center. While there are a number of different methods to discretize the spatial dimension for both one- and two-dimensional geometries, this work focuses on the second order diamond difference method [LM93]. The first step in this process is to break down the streaming operator (\(\bsOmega \cdot \nabla \psi_{m, g}\)) for rectangular and spherical coordinates. From there, the diamond difference will be used to solve for the flux at specific cell edges and average the neutron flux at the cell center through
where \(i\) and \(j\) represent the \(x\) and \(y\) spatial cells, respectively with \(i \in [1, \ldots, I]\) and \(j \in [1, \ldots, J]\) and the \(\psi_{i \pm 1/2}\) flux values representing the spatial cell edges. If the incoming flux is known at \(i - 1/2\), the average flux for each spatial cell can be solved through the spatial sweep. There is a different spatial discretization method that will be briefly used in a later verification chapter, which is a first order piece-wise constant upwind scheme known as the step method. Unless noted otherwise, the diamond difference is the spatial discretization scheme employed in this analysis.
It should be noted that for clarity, fixed source problems will be discussed in this section and both the time-dependent term and the \(g\) subscripts will be dropped from Equation (1). This alters the neutron transport equation for a steady-state, mono-energetic problem as
which will be the basis for the one- and two-dimensional spatial discretizations of the streaming operator.
2.3.1. One-Dimensional Slab#
The simplest neutron transport is a one-dimensional slab, in which the streaming operator is represented as
where \(x\) is the spatial dimension and \(\mu\) is the angular dimension. The partial derivative with respect to \(x\) can be approximated with the fundamental theorem of calculus by taking
where the denominator can be simplified as \(\Delta x_{i} = x_{i+1/2} - x_{i-1/2}\). Substituting this back into the neutron transport equation, the equation becomes
when solving for a specific spatial cell \(i\) and angle \(\mu_{m}\).
The diamond difference method, defined in Equation (1) can replace the angular flux in the removal term in which the equation becomes
The sign of the angular dimension \(\mu\) determines the angular flux to solve for, in which the incoming boundary flux is known at \(\psi_{1/2,m}\) when \(\mu_{m} > 0\) and likewise at \(\psi_{I+1/2,m}\) when \(\mu_{m} < 0\). In the instance where the angle \(\mu_{m} > 0\), the angular flux at \(\psi_{i-1/2,m}\) is known and the angular flux at \(\psi_{i+1/2,m}\) is calculated as
where the \(\ell\) and \(\ell+1\) are added to show an iterative approach to solving this equation with the scalar flux being lagged by an iteration. The scalar flux is updated according to Equation (2).
2.3.2. One-Dimensional Sphere#
The spatial discretization of the one-dimensional spherical geometry is more involved, as the streaming operator becomes
where \(x\) is the spatial dimension and \(\mu\) is the angular dimension. Following [LM93], angular differencing coefficients \(\alpha\) are added to the transport equation to become
where \(w_{m}\) is the discrete ordinates weight for angle \(m\). This introduces the half angle \(m \pm 1/2\) term where
is the diamond difference method in the angular dimension. The half angle is calculated as
and the starting and ending directions are \(\mu_{1/2} = -1\) and \(\mu_{M+1/2} = 1\). For the streaming term to vanish in a uniform flux problem, the angular differencing coefficient must have the relationship
and the property \(\alpha_{1/2} = \alpha_{M+1/2} = 0\) to satisfy the neutron balance condition in spherical geometry. Substituting Equation (8) into Equation (7) for the \(\psi_{m+1/2}(x)\) term, the transport equation becomes
where \(\psi_{m}(x)\) can be solved when \(\psi_{m-1/2}(x)\) is known.
Before proceeding with the spatial sweep, the partial derivative with respect to \(x\) in Equation (11) must be handled appropriately in curvilinear coordinates. For this, the volume \(V_{i}\) of one spherical shell is
where \(i \pm 1/2\) are the cell edges and the surface area at the cell boundary is \(A_{i+1/2} = 4 \, \pi \, x^2_{i+1/2}\). Integrating Equation (11) by Equation (12) and dividing by the volume, the spherical neutron transport equation becomes
for a given angle \(m\) and spatial cell \(i\).
With Equation (13), the half angular flux is required to calculate the flux at a given angle \(m\) and cell \(i\). To do this, the angular derivative in Equation (6) can be used with the starting direction \(\mu_{1/2} = -1\) as the inward most direction to cancel out the angular derivative. The neutron transport equation for spherical geometry would simplify to
which can be rearranged to solve for the starting angular flux \(\psi_{1/2,i}\) as
with \(\Delta x_{i} = x_{i+1/2} - x_{i-1/2}\).
After the flux for the half angle \(m = 1/2\) is calculated, Equation (13) can be rearranged to solve for the cell flux center as
for \(\mu < 0\) with
where flux at \(\psi_{m,i+1/2}\) is known. The \(\psi_{m,i-1/2}\) is updated using the diamond difference method as
and the angular flux half angle \(\psi_{m+1/2,i}\) can be updated using Equation (8). If sweeping from the center of the sphere to the edge, \(\mu > 0\), the \(a_{m,i}\) and \(d_{m,i}\) equations in Equation (16) are updated as
where the angular flux \(\psi_{m, i+1/2}\) is updated using the diamond difference method.
The center of the sphere can approximated as
and treated as a reflective slab boundary condition. This however introduces a nonphysical dip in the spatial flux as pointed out in [LM93] and it would be more appropriate to approximate the flux as
where the starting flux is calculated using Equation (14) at the origin for all angles. This can be combined with a weighted diamond difference scheme as presented in [MM84] to update the half angle as
instead of using Equation (8). The weight \(\tau\) is
where the \(\mu_{m \pm 1/2}\) are updated according to Equation (9).
All the equations are presented to solve for the one-dimensional sphere problem using discrete ordinates and the diamond difference spatial discretization. Assuming the angular dimension is traversed from \(\mu = -1\) to \(\mu = 1\), the initial conditions are \(\mu_{1/2} = -1\) and \(\alpha_{1/2} = 0\), which can solve for the half angle \(\psi_{1/2,i}\) using Equation (14). The angular flux at the spatial cell centers using Equation (16), where the spatial cell edge is updated using Equation (17) and the angle edge is updated using Equation (20). The angular differencing coefficient and half angle term are updated using Equations (10) and (9). When the spatial sweep is from the spherical center to the edge \(\mu > 0\), the half angle Equation (19) is used for the initialization and the \(a_{m,i}\) and \(d_{m,i}\) terms in Equation (15) are updated with Equation (18).
2.3.3. Two-Dimensional Slab#
The streaming term in two dimensions would result in
where \(\mu\) and \(\eta\) are the angular dimensions and \(x\) and \(y\) are the spatial dimensions. As with the one-dimensional slab geometry, the fundamental theorem of calculus can be used to discretize the partial derivative with respect to \(y\) as
with the derivative for \(x\) shown in Equation (5) and \(\Delta y_{j} = y_{j+1/2} - y_{j-1/2}\) equaling the cell width in the \(y\) direction. Inserting these terms back into the neutron transport equation, the resulting equation is
when solving for a specific spatial cell \((i, j)\) and angle \(\mu_{m}, \eta_{m}\).
For the specific discrete ordinate, the same principle used for the known flux in the one-dimensional slab problem can be expanded to the two-dimensional slab. Assuming that the discrete ordinate results in \(\mu_{m} > 0\) and \(\eta_{m} > 0\), the fluxes at the boundaries \(\psi_{i-1/2,j,m}\) and \(\psi_{i,j-1/2,m}\) will be known. The diamond difference method in Equations (2) and (3) can be rearranged to solve for the opposite sides, namely \(\psi_{i+1/2,j,m}\) and \(\psi_{i,j+1/2,m}\). These terms are reinserted into the transport equation in which only the angular flux at the cell center \(\psi_{i,j,m}\) is unknown. Combining like terms, the two-dimensional transport sweep equation is
which can be solved in the same iterative scheme as the one-dimensional slab case. The angular flux at the cell edges, \(\psi_{i+1/2,j,m}\) and \(\psi_{i,j+1/2,m}\), are updated using the diamond difference Equations (2) and (3) before proceeding to the next spatial cell.