Table of Contents
This page is under construction. If it seems incomplete, it is, but feel free to enjoy what is already here. This page has an accompanying toy nodal expansion method code, which lives here.
Transport Equation
The neutron transport equation defines the balance of neutron creation and neutron removal processes from a volume $V$. For a system with no external source, the time rate of change of the flux is given by $$ \begin{align} \frac{1}{v}\frac{\partial\phi}{\partial t}&= \int\limits_{4\pi}\int\limits_0^\infty \Sigma_s(\vec r,E’\rightarrow E,\hat\Omega’\rightarrow\hat\Omega,t) \phi(\vec r,E’,\hat\Omega’,t)dE’d\hat\Omega’ \\ &+\frac{\chi(E)}{4\pi}\int\limits_0^\infty \nu(E’)\Sigma_f(\vec r,E’,t)\phi(\vec r,E’,\hat\Omega,t)dE’ \\ &+\frac{1}{4\pi}\sum\limits_f\chi_f(E)\lambda_fC_f(\vec r,t) \\ &-\hat\Omega\cdot\nabla\phi(\vec r,E,\hat\Omega,t) \\ &-\Sigma_t(\vec r,E,t)\phi(\vec r,E,\hat\Omega,t) \end{align} $$ where
- $v$ is the neutron velocity
- $\phi$ is the neutron flux, and is a function of position, energy, angle, and time
- $\Sigma_s$ is the scattering cross section, and is a function of position, energy, angle, and time
- $\chi(E)$ is the probability a neutron produced from fission was produced at energy $E$
- $nu$ is the average number of neutrons produced in a fission event, and is a function of incident neutron energy
- $\Sigma_f$ is the fission cross section, and is a function of position, energy, and time
- $\chi_f(E)$ is the probability a decay neutron from fission product $f$ is emitted at energy $E$
- $\lambda_f$ is the decay constant of fission product $f$
- $C_f$ is the concentration of fission product $f$, and is a function of position and time
- $\Sigma_t$ is the total cross section, and is a function of position, neutron energy, and time
On the right-hand side of the equation, there are positive terms that represent neutron gains in the system, and negative terms that represent neutron losses in the system. At this point, let us consider the flux of neutrons with a fixed energy $E$ and direction $\hat\Omega$ at the point $\vec r$ in space. The first term on the right-hand side of the transport equation is the in-scattering source, where neutrons at a different angle and energy scatter into our energy $E$ and angle $\hat\Omega$. The next term describes the neutrons produced at our position, angle, and energy from fission. The final gain term is the production of neutrons from the decay of radioactive fission products. The first loss term is the leakage of neutrons away from our position $\vec r$. The other loss term simply represents neutrons being removed from our position, angle, or energy from suffering an interaction in matter.
Discretization of the Transport Equation
The full time-dependent neutron transport equation is very difficult to solve. There are many different approximations to the transport equation that are available which allow it to be calculated. One such approximation is to assume that the angular distribution does not strongly effect the solution. To do this, we integrate the transport equation over all angles to get $$ \begin{align} \frac{1}{v}\frac{\partial\phi}{\partial t}&= \int\limits_0^\infty\Sigma_s(E’\rightarrow E)\phi(\vec r,E’,t) \\ &+\frac{\chi(E)}{4\pi}\int\limits_0^\infty \nu(E’)\Sigma_f(\vec r,E’,t)\phi(\vec r,E’,t)dE’ \\ &+\frac{1}{4\pi}\sum\limits_f\chi_f(E)\lambda_fC_f(\vec r,t) \\ &-\int\limits_{4\pi} \hat\Omega\cdot\nabla\phi(\vec r,E,\hat\Omega,t)d\hat\Omega \\ &-\Sigma_t(\vec r,E,t)\phi(\vec r,E,t) \end{align} $$ We can then use the divergence theorem on the first loss term, which yields the neutron continuity equation: $$ \begin{align} \frac{1}{v}\frac{\partial\phi}{\partial t}&= \int\limits_0^\infty\Sigma_s(E’\rightarrow E)\phi(\vec r,E’,t) \\ &+\frac{\chi(E)}{4\pi}\int\limits_0^\infty \nu(E’)\Sigma_f(\vec r,E’,t)\phi(\vec r,E’,t)dE’ \\ &+\frac{1}{4\pi}\sum\limits_f\chi_f(E)\lambda_fC_f(\vec r,t) \\ &-\nabla\cdot\vec J(\vec r,E,t) \\ &-\Sigma_t(\vec r,E,t)\phi(\vec r,E,t) \end{align} $$ where $\vec J$ is the neutron current.
Next, we discretize the equation in energy. We say that our energy spectrum is divided up into $G$ groups, where the properties within each group $g$ are approximated to be constant. This allows us to rewrite the continuity equation as $$ \frac{1}{v}\frac{\partial\phi_g}{\partial t} +\nabla\cdot\vec J_g(\vec r,t) +(\Sigma_{t,g}(\vec r,t)-\Sigma_{s,g}(\vec r,t))\phi_g(\vec r,t) =Q_g(\vec r,t) $$ where the gain terms have been combined into a single source term $Q$.
At this point, we will choose to solve a steady-state problem, meaning that the time dependence of the terms can be removed. In order to solve this, we also need to discretize our problem in space. Dividing space into cubes turns out to be convenient. Each cube in space is approximated to have constant properties. We call one of these cubes a node, and denote a specific node as $k$. This leads us to the following form of the diffusion equation: $$ \nabla\cdot\vec J^k_g+\Sigma^k_{r,g}\phi^k_g=Q^k_g $$ where we use the removal cross section, defined as $$ \Sigma^k_{r,g}=\Sigma^k_{a,g} +\sum\limits^G_{g’\neq g}\Sigma^k_{s,g\rightarrow g’} $$ and $$ Q^k_g(\vec{r})=\frac{\chi_g}{k_{eff}} \sum\limits_{g’=1}^G\nu\Sigma_{f,g’}^k\phi^k_{g’}(\vec{r})+ \sum\limits_{g’\neq g}^G\Sigma^k_{s,g’\rightarrow g}\phi^k_{g’}(\vec{r}) $$ This is the starting point of the derivation in Okumura, 1998, which is roughly translated in the following section.
The Diffusion Approximation
We now introduce another approximation by assuming that the relationship between the flux and the current follow Fick’s Law. This can be expressed as $$ \vec J_g(\vec r,t)=-D_g(\vec r)\nabla\phi_g(\vec r,t) $$ We will use Fick’s Law in the following sections to eliminate an unknown (the current) in a system of equations, a common result in an introductory reactor physics course.
As an aside, plugging this into the continuity equation gives the diffusion equation: $$ \frac{1}{v}\frac{\partial\phi}{\partial t}- \nabla\cdot D_g(\vec r)\nabla\phi_g(\vec r,t) +\Sigma_{a,g}(\vec r,t)\phi_g(\vec r,t)=Q_g(\vec r,t) $$
Nodal Balance
We define the center of our node $k$ to be $(0,0,0)$. Therefore, the node spans the space $-\frac{\Delta x}{2}<x<\frac{\Delta x}{2}$, $-\frac{\Delta y}{2}<y<\frac{\Delta y}{2}$, $-\frac{\Delta z}{2}<z<\frac{\Delta z}{2}$. “+” and “-” are used to denote the face of the node with a positive or negative component, respectively. Now, we expand the divergence of the current in the diffusion equation, integrate over the node, and divide by the node volume. This results in the node balance equation, given by $$ \frac{1}{\Delta_x^k}\left(J^k_{g,x+}-J^k_{g,x-}\right)+ \frac{1}{\Delta_y^k}\left(J^k_{g,y+}-J^k_{g,y-}\right)+ \frac{1}{\Delta_z^k}\left(J^k_{g,z+}-J^k_{g,z-}\right)+ \Sigma^k_{rg}\phi^k_{g,0}=Q^k_{g,0} $$ where we defined $$ \phi^k_{g,0}=\frac{1}{\Delta^k_x\Delta^k_y\Delta^k_z} \int\limits^{\Delta^k_x/2}_{-\Delta^k_x/2} \int\limits^{\Delta^k_y/2}_{-\Delta^k_y/2} \int\limits^{\Delta^k_z/2}_{-\Delta^k_z/2} \phi^k_g(x,y,z)dxdydz $$ and $$ J^k_{g,x\pm}=\frac{1}{\Delta^k_y\Delta^k_z} \int\limits^{\Delta^k_y/2}_{-\Delta^k_y/2} \int\limits^{\Delta^k_z/2}_{-\Delta^k_z/2} J^k_g\left(\pm\frac{\Delta^k_x}{2},y,z\right)dydz $$ and $$ \begin{align} Q^k_{g,0}&=\frac{1}{\Delta^k_x\Delta^k_y\Delta^k_z} \int\limits^{\Delta^k_x/2}_{-\Delta^k_x/2} \int\limits^{\Delta^k_y/2}_{-\Delta^k_y/2} \int\limits^{\Delta^k_z/2}_{-\Delta^k_z/2} Q^k_g(x,y,z)dxdydz \\ &=\frac{1}{\Delta^k_x\Delta^k_y\Delta^k_z} \int\limits^{\Delta^k_x/2}_{-\Delta^k_x/2} \int\limits^{\Delta^k_y/2}_{-\Delta^k_y/2} \int\limits^{\Delta^k_z/2}_{-\Delta^k_z/2} \left(\frac{\chi_g}{k_{eff}} \sum\limits_{g’=1}^G\nu\Sigma_{fg’}^k\phi^k_{g’}(x,y,z)+ \sum\limits_{g’\neq g}^G\Sigma^k_{sg’\rightarrow g}\phi^k_{g’}(x,y,z) \right)dxdydz \\ &=\frac{\chi_g}{k_{eff}} \sum\limits_{g’=1}^G\nu\Sigma_{fg’}^k\phi^k_{g’,0}+ \sum\limits_{g’\neq g}^G\Sigma^k_{sg’\rightarrow g}\phi^k_{g’,0} \end{align} $$ The $y$ and $z$ currents are defined in the same way as the $x$ current. In order to reduce the number of variables, we would like to reduce the number of unknowns. One approximation often made to neutron transport is the application of Fick’s Law, which relates the current to the flux as $$ \vec J^k_g=-D^k_g\nabla\phi^k_g $$ Plugging this into the definition of the $x$ current yields $$ J^k_{g,x\pm}=\frac{1}{\Delta^k_y\Delta^k_z} \int\limits^{\Delta^k_y/2}_{-\Delta^k_y/2} \int\limits^{\Delta^k_z/2}_{-\Delta^k_z/2} -D^k_{g,x}\frac{\partial}{\partial x}\phi^k_g(x,y,z) \biggr|_{x=\pm\Delta^k_x/2}dydz $$
The node balance equation is expressed in terms of the volume-averaged neutron flux within a node, denoted as $\phi^k_{g,0}$, and the surface-averaged net neutron currents across the six faces of the node. To proceed, we introduce the concepts of surface-averaged neutron flux and surface-averaged partial neutron currents. Within diffusion theory, the partial neutron currents are typically defined using the P1 approximation of the angular neutron flux, as shown below. $$ \begin{align} J^{out}_g(\vec{r})&=\frac{\phi_g(\vec{r})}{4}- \frac{D_g}{2}\hat{e}_s\cdot\nabla\phi(\vec{r}) \\ J^{in}_g(\vec{r})&=\frac{\phi_g(\vec{r})}{4}+ \frac{D_g}{2}\hat{e}_s\cdot\nabla\phi(\vec{r}) \end{align} $$ Here, $\hat{e}_s$ represents the unit outward normal vector to the surface under consideration. The outgoing partial neutron current in the direction of this outward normal vector is denoted as $J^{out}$, while the incoming partial neutron current, flowing in the opposite direction, is denoted as $J^{in}$. The partial current leaving the x+ surface is therefore $$ J^{out,k}_{g,x+}=\frac{1}{4}\phi^k_{g,x+}+\frac{1}{2} \left[ \frac{1}{\Delta^k_y\Delta^k_z} \int\limits^{\Delta^k_y/2}_{-\Delta^k_y/2} \int\limits^{\Delta^k_z/2}_{-\Delta^k_z/2} -D^k_{g,x}\frac{\partial}{\partial x}\phi^k_g(x,y,z) \biggr|_{x=\Delta^k_x/2}dydz \right] $$ where $\phi^k_{g,x+}$ is the surface-averaged neutron flux on the $k+$ interface of the positive $x$ surface of the node. $$ \phi_{g,x+}^{k}=\frac{1}{\Delta_y^k\Delta_z^k} \int\limits_{-\Delta_y^k/2}^{\Delta_y^k/2} \int\limits_{-\Delta_z^k/2}^{\Delta_z^k/2} \phi^{k}_g\left(\frac{\Delta_x^k}{2},y,z\right)dydz $$
The term in the square brackets is the definition of $J^k_{g,x+}$ from above. Substituting this in gives that $$ J^{out,k}_{g,x+}=\frac{1}{4}\phi^k_{g,x+}+\frac{1}{2}J^k_{g,x+} \ $$ Similarly, since the only difference in the definition of the incoming current from the P1 approximation is the sign flip, we can show that $$ J^{in,k}_{g,x+}=\frac{1}{4}\phi^k_{g,x+}-\frac{1}{2}J^k_{g,x+} $$ Adding and subtracting these two equations gives expressions for the current and flux at the $x+$ surface of the node: $$ \begin{align} J^k_{g,x+}&=J^{out,k}_{g,x+}-J^{in,k}_{g,x+} \\ \phi^k_{g,x+}&=2(J^{out,k}_{g,x+}+J^{in,k}_{g,x+}) \end{align} $$ Mirroring these calculations for the $x-$ surface yields $$ \begin{align} J^k_{g,x-}&=J^{in,k}_{g,x-}-J^{out,k}_{g,x-} \\ \phi^k_{g,x-}&=2(J^{out,k}_{g,x-}+J^{in,k}_{g,x-}) \end{align} $$ This allows us to rewrite the node balance equation using the incoming and outgoing currents. $$ \begin{align} \frac{1}{\Delta^k_x} (J^{out,k}_{g,x+}-J^{in,k}_{g,x+}-J^{in,k}_{g,x-}+J^{out,k}_{g,x-}) &+\frac{1}{\Delta^k_y} (J^{out,k}_{g,y+}-J^{in,k}_{g,y+}-J^{in,k}_{g,y-}+J^{out,k}_{g,y-}) \nonumber \\ &+\frac{1}{\Delta^k_z} (J^{out,k}_{g,z+}-J^{in,k}_{g,z+}-J^{in,k}_{g,z-}+J^{out,k}_{g,z-})+ \Sigma^k_{r,g}\phi^k_{g,0}=Q^k_{g,0} \end{align} $$
We do this because it allows us to split up the calculation. I will show in a later section that we can solve for the outgoing current on each surface of a node from the values of the incoming currents. This allows a two-step iteration strategy. First, we calculate the outgoing current for each node from its incoming current. The outgoing current for each node is used to update the incoming currents of the adjacent nodes. We can iterate this until the system of nodes is balanced to converge on a solution.
Transverse Integration
In order to reduce the solution complexity, we derive a one dimensional diffusion equation by integrating over the transverse directions. This is performed for each direction, yielding three equations toward the solution. The steps to derive the diffusion equation in x are roughly outlined in this section. The same steps can be equivalently followed for integration in the y and z directions as well.
Earlier, we derived a special form of the diffusion equation. $$ \nabla\cdot\vec J^k_g+\Sigma^k_{r,g}\phi^k_g=Q^k_g $$ We can expand the divergence of the current and integrate with respect to y and z within the node. Dividing by the area $\Delta_y^k\Delta_z^k$ then yields the following: $$ \begin{align} &\frac{1}{\Delta^k_y\Delta^k_z} \int\limits_{-\Delta_y^k/2}^{\Delta_y^k/2} \int\limits_{-\Delta_z^k/2}^{\Delta_z^k/2} \frac{\partial}{\partial x}J^k_{g,x}(x,y,z)dydz +\frac{1}{\Delta^k_y\Delta^k_z} \int\limits_{-\Delta_y^k/2}^{\Delta_y^k/2} \int\limits_{-\Delta_z^k/2}^{\Delta_z^k/2} \frac{\partial}{\partial y}J^k_{g,y}(x,y,z)dydz \nonumber \\ +&\frac{1}{\Delta^k_y\Delta^k_z} \int\limits_{-\Delta_y^k/2}^{\Delta_y^k/2} \int\limits_{-\Delta_z^k/2}^{\Delta_z^k/2} \frac{\partial}{\partial z}J^k_{g,z}(x,y,z)dydz +\frac{1}{\Delta^k_y\Delta^k_z} \int\limits_{-\Delta_y^k/2}^{\Delta_y^k/2} \int\limits_{-\Delta_z^k/2}^{\Delta_z^k/2} \Sigma^k_{r,g}\phi^k_g(x,y,z)dydz \nonumber \\ =&\frac{1}{\Delta^k_y\Delta^k_z} \int\limits_{-\Delta_y^k/2}^{\Delta_y^k/2} \int\limits_{-\Delta_z^k/2}^{\Delta_z^k/2} Q^k_g(x,y,z)dydz \end{align} $$ which can be simplified to: $$ \begin{equation} \frac{d}{dx}\bar J^k_{g,x}(x)+\Sigma^k_{r,g}\bar\phi^k_{g,x}(x) =\bar Q^k_{g,x}(x)-\frac{1}{\Delta^k_y}\bar L^k_{g,y}(x) -\frac{1}{\Delta^k_z}\bar L^k_{g,z}(x) \end{equation} $$ where $$ \begin{align} \bar\phi^{k}_{g,x}(x)&=\frac{1}{\Delta^k_y\Delta^k_z} \int\limits_{-\Delta_y^k/2}^{\Delta_y^k/2} \int\limits_{-\Delta_z^k/2}^{\Delta_z^k/2} \phi^k_g(x,y,z)dydz \nonumber \\ \bar J^{k}_{g,x}(x)&=\frac{1}{\Delta^k_y\Delta^k_z} \int\limits_{-\Delta_y^k/2}^{\Delta_y^k/2} \int\limits_{-\Delta_z^k/2}^{\Delta_z^k/2} J^k_g(x,y,z)dydz \nonumber \\ \bar Q^k_{g,x}(x)&=\frac{1}{\Delta^k_y\Delta^k_z} \int\limits_{-\Delta_y^k/2}^{\Delta_y^k/2} \int\limits_{-\Delta_z^k/2}^{\Delta_z^k/2} Q^k_g(x,y,z)dydz \nonumber \\ &=\frac{\chi^k_g}{k} \sum\limits_{g’=1}^{G}\nu\Sigma^{k}_{f,g’}\bar\phi^k_{g’}(x)+ \sum\limits_{g’\neq g}^{G}\Sigma^{k}_{s,g’\rightarrow g} \bar\phi^k_{g’}(x) \\ \bar L^{k}_{g,y}(x)&=\frac{1}{\Delta^k_z} \int\limits_{-\Delta_z^k/2}^{\Delta_z^k/2}\left[ -D^k_{g,y}\frac{\partial\phi^k_g(x,y,z)}{\partial y} \right]^{\Delta^k_y/2}_{-\Delta^k_y/2}dz \nonumber \\ \bar L^{k}_{g,z}(x)&=\frac{1}{\Delta^k_y} \int\limits_{-\Delta_{y}^{k}/2}^{\Delta_{y}^{k}/2}\left[ -D^k_{g,z}\frac{\partial\phi^{k}_{g}(x,y,z)}{\partial z} \right]^{\Delta^k_z/2}_{-\Delta^k_z/2}dy \end{align} $$ are the parameters averaged over the area of the node.
Similar equations can be written for the y and z directions. This is the transverse-integrated one-dimensional neutron balance equation. The equation is a one-dimensional diffusion equation with transverse leakage terms.
We can now write Fick’s Law in the x-direction averaged over a node. $$ \begin{equation} \bar J^k_{g,x}(x)=-D^k_{g,x}\frac{d}{dx}\bar\phi^k_{g,x}(x) \end{equation} $$ and similarly for the y and z directions. Substituting these into the 1D diffusion equation with transverse leakage terms results in the following one-dimensional diffusion equation. $$ \begin{equation} -D^k_{g,x}\frac{d^2}{dx^2}\bar\phi^k_{g,x}(x)+ \Sigma^k_{r,g}\bar\phi^k_{g,x}(x) =\bar Q^k_{g,x}(x)-\frac{1}{\Delta^k_y}\bar L^k_{g,y}(x) -\frac{1}{\Delta^k_z}\bar L^k_{g,z}(x) \end{equation} $$