Solving the Time-Independent Diffusion Equation

Quick foreword: This is part two of a guide that serves as both a refresher of a class I very much enjoyed taking by Professor Massimiliano Fratoni and as a resource on something that I think is a bit harder to understand from the ground up. This part will be a bit more math heavy and will reference a lot more multivariable calculus and differential equations than the previous part. I'll still try to keep it somewhat palatable for anyone who's read the last part.

In the previous part of this series, we derived the Time-Independent Diffusion Equation, but we didn't fully solve it out. We'll first do that here in this section so that we can gain some insight when we're looking at reactors of different shapes.

As a reminder, the Time-Independent Diffusion Equation (TIDE) can be written as:

Time-Independent Diffusion Equation

$$0 = S - \Sigma_{abs}\,\phi + D\,\nabla^2\phi$$

Where:

$S$the source term (where our neutrons come from)

$\Sigma_{abs}$the macroscopic neutron absorption cross section

$\phi$our neutron flux

$D$our diffusion coefficient

(If you need a review on what any of these things are, please visit the last part of the section here.)

Part 1: Boundary Conditions

In order for our TIDE to be physically accurate, we have some basic conditions that our flux must satisfy. To start with, our $\phi$ must be positive, real and finite. These conditions are trivially required as our system would be unphysical without these properties (if a reactor produced an infinite amount of flux, it would be pretty catastrophic and would violate the laws of physics).

The other set of primary conditions concerns interfaces, an interface is any area where once crossed, requires a different set of parameters to continue describing physical behavior. For example, if we are considering how neutrons flow through a shielding wall, there exists a wall-air interface. At these interfaces, we need to make sure that the predicted flux by both sides are the same. In the wall-air example, this means our model of neutron flux through a wall has to predict the same value at the interface of neutron flux through the air.

We can generalize this by saying for any interfaces $A$ and $B$ at position $L$, the following holds:

Flux is continuous across an interface

$$\phi_A(L) = \phi_B(L)$$

This rule also holds for currents at interfaces. If you recall from the last part, the current $J(x) = -D\,\nabla\phi$.

This is a result of the limit definition of a derivative, we know that the derivative must equal the following limit:

$$f'(x) = \lim_{h \to 0} \frac{f(x-h) - f(x)}{h}$$

If our derivatives from both sides are not equal at the interface, it means that the flux reaches some discontinuity at that point. This obviously is not physically possible, so we must constrain the boundary as follows:

Current is continuous across an interface

$$J_A(L) = J_B(L)$$

Addendum 1: Fick's Law Correction

Readers familiar with Fick's law may take issue with the behavior of this law at the boundaries of the geometries that we choose. Fick's law has four necessary assumptions it makes to work (for this article I will not be explaining them, but there's plenty of reading on it if needed):

The flux must be almost the same regardless of whatever angle we measure it at. There can be a slight variance in a specific direction which we define as the current.
Fluxes shouldn't vary too much within short distance, that is we need gradual changes in flux over distances
We assume that currents immediately change when the flux gradient changes, in reality there is some inertia there but in order to apply Fick's law we must assume it is negligible
You must have weak absorption, this means you cannot have any sort of absorber which absorbs at a rate where it would create neutron flux anisotropy (this is a special case of 1)

Sharp readers may notice that these conditions aren't met at a boundary as boundaries create strongly directional neutron flux conditions. In order to still use Fick's Law, we have some extrapolation distance $d$ beyond the surface where the flux vanishes at. We will not be using it in this case as it is negligible for our uses but for reference it is good to know that

Extrapolation distance

$$d = 2.13\,D$$

Part 2: Solving the Differential Equations

So now we have some general conditions established by our boundary conditions, we can now solve the differential equation. In order to solve our equations, we need to first define the geometry of the source we aim to describe. This geometry will affect the specific solution we get out of the diffusion equation.

In this article, we will discuss solving for two separate geometries: the infinite planar source and the point source. In further sections, we will use similar techniques to explore reactors of various geometries. Keep in mind that the planar source will be a bit simpler to derive and the point source may add a bit of complexity. Do not worry if you are unable to tackle the point source, we will cover it in further depth in future sections, it is merely included as a demonstration here.

Infinite Planar Source

An infinite planar source is essentially an infinitely large, infinitely thin sheet that has a constant flux across the $y$ and $z$ directions. For this system, we are going to center our planar source on the $x$-axis by placing it at $x = 0$, with there being no other sources at any other $x$.

FIG. 1 — INFINITE PLANAR SOURCE. A thin sheet at x = 0 emits neutrons isotropically into both half-spaces. Because nothing varies along y or z, the problem collapses to one dimension. The flux peaks at the source and falls off exponentially with distance, symmetric about x = 0.

Now, if we look at the area outside of our source (every $x$ coordinate that is not zero), we know that there is no source term for our differential equation. As a result, we can write the diffusion equation at all nonzero points to be:

Planar diffusion equation, valid for x ≠ 0

$$\frac{d^2 \phi}{dx^2} - \frac{1}{L^2}\,\phi = 0$$

(We only differentiate with respect to $x$ as the gradient in the $y$ and $z$ directions is definitionally zero due to our constant flux in those directions.)

If you have taken a differential equations class before, you'll recognize this is a second order differential equation¹×See: tutorial.math.lamar.edu/classes/de/introsecondorder.aspx. As such, it has the general solution:

General solution

$$\phi(x) = A\,e^{-x/L} + C\,e^{x/L}$$

You can see now that we have two unknown constants $A$ and $C$. From looking at the equation, we can immediately point the $C\,e^{x/L}$ part out as a red flag. We can see that flux increases with distance from the source, which would not match our physical model for how these sources actually work. In addition, this allows for essentially infinite flux at arbitrarily large distances meaning it violates our original conditions of not having infinite flux. From this we can conclude that in order for the solution to be valid $C$ must equal $0$.

Now we are left with:

$$\phi(x) = A\,e^{-x/L}$$

In order to find our unknown constant $A$, we can use some reasoning about the nature of our source. We know that the planar source shoots out neutrons isotropically meaning it equally shoots out neutrons in both directions. This means that the total flux of our source $S$ will be equal to the sum of the current going into the positive and negative directions. Due to the fact that our source is even on both sides, we know that the total flux $S$ is just two times our current at the source. Therefore:

$$J(0) = \frac{S}{2}$$

Now, let's go back to our definition of current. We know the current $J$ to be:

$$J(0) = -D\,\nabla\phi(0)$$

Let's plug in the $\phi$ we got from above:

$$J(0) = -D\,\frac{d}{dx}\Big|_{x=0}\!\left(A\,e^{-x/L}\right)$$

If we take our derivative and evaluate at $0$ we get the following:

$$J(0) = \frac{DA}{L}$$

Now we have two forms for $J(0)$ which allow us to solve for our unknown $A$ in terms of our known constants $S$, $D$ and $L$:

$$\frac{S}{2} = \frac{DA}{L} \quad\Longrightarrow\quad A = \frac{SL}{2D}$$

If we plug that back in to our original equation for $\phi$, we get our final form of:

Infinite planar source — neutron flux

$$\phi(x) = \frac{SL}{2D}\,e^{-|x|/L}$$

We include the absolute value for the $x$ as the flux only changes with distance from the source, it doesn't matter if the distance is positive or negative.

Point Source

FIG. 2 — POINT SOURCE. A single point emits neutrons isotropically, so the flux depends only on radial distance r. Spreading the same neutrons over ever-larger spherical shells (the 1/r factor) together with absorption (the e^(−r/L) factor) makes the flux fall off sharply — and diverge as r → 0.

For the point source, we're going to imagine a single point in space that is an isotropic source of neutrons. Because of its isotropic nature, this is much easier to model in a spherical coordinate system²×See: tutorial.math.lamar.edu/classes/calciii/sphericalcoords.aspx. We can rewrite our diffusion equation spherically as follows:

Diffusion equation in spherical coordinates

$$\frac{1}{r^2}\frac{d}{dr}\!\left(r^2 \frac{d\phi}{dr}\right) - \frac{1}{L^2}\,\phi = 0$$

Similar to the planar source, this equation holds whenever there isn't a source, which means any point that isn't the point involved in our singular point source.

Now we're going to define a dummy variable $w$ in order to make solving this differential equation easier by setting $w = r\phi$. Using this we get the following form:

$$\frac{d^2 w}{dr^2} - \frac{1}{L^2}\,w = 0$$

We now have a differential equation we can solve exactly the same as we did last time, as the following:

$$w = A\,e^{-r/L} + C\,e^{r/L}$$

We already know from the planar source that $C$ must equal zero, so we are left with:

$$w = A\,e^{-r/L} \quad\Longrightarrow\quad \phi = \frac{A}{r}\,e^{-r/L}$$

We can find $A$ using our current again: we need to split the total source neutrons evenly across every direction to find our current $J$. Doing that we get that:

$$J(0) = \lim_{r \to 0} \frac{S}{4\pi r^2}$$

Now conducting the exact same steps as our previous section by plugging $\phi$ into our current formula and solving for $A$, we get:

$$A = \frac{S}{4\pi D}$$

Plugging that back into our $\phi$ just as before gets us:

Point source — neutron flux

$$\phi(r) = \frac{S\,e^{-r/L}}{4\pi D r}$$

This is our final solution for the point source.

Both solutions share the same skeleton: an exponential decay set by the diffusion length³×L is the diffusion length — the characteristic distance over which the flux falls to 1/e (about 37%) of its value. Larger L means neutrons travel further before being absorbed. $L$. The geometry only changes the prefactor: a constant for the plane, a $1/r$ falloff for the point. Drag the control below to feel how $L$ sets the reach of the flux.

Diffusion length explorer

Drag left/right across the plot — or use the slider — to change the diffusion length L. Switch the source geometry to compare how the flux reaches into space. The dashed marker sits at one diffusion length, where the flux has dropped to 1/e.

Diffusion length L1.0 cm

In the next part of this series, we will explore diffusion in nuclear reactors and how to model it starting with one-group reactor theory. The concepts in this part and the previous part will be very useful for understanding this material.

Section Notes

¹ Second-order ODEs: see tutorial.math.lamar.edu/classes/de/introsecondorder.aspx

² Spherical coordinates: see tutorial.math.lamar.edu/classes/calciii/sphericalcoords.aspx

³ L is the diffusion length — the characteristic distance over which the flux falls to 1/e of its value. Larger L means neutrons travel further before being absorbed.

Back · Part 1 Neutron Diffusion From the Ground Up ← The neutron, cross sections, and deriving the diffusion equation

AI Usage: Claude was used to help code the diagrams and the explorer from mockups I made; the content is all mine!

← back to essays