What is Diffusion Equation – Definition

The diffusion equation is mostly solved in media with high densities such as neutron moderators (H₂O, D₂O or graphite). The problem is usually bounded by air. The mean free path of neutron in air is much larger than in the moderator, so that it is possible to treat it as a vacuum in neutron flux distribution calculations. The vacuum boundary condition supposes that no neutrons are entering a surface.

If we consider that no neutrons are reflected from the vacuum back to the volume, the following condition can be derived from the Fick’s law:

Where d ≈ ⅔ λ_tr is known as the extrapolated length. For homogeneous, weakly absorbing media, an exact solution of the mono-energetic transport equation in this case yields d ≈ 0.7104 λ_tr. The geometric interpretation of the previous equation is that the relative neutron flux near the boundary has a slope of -1/d, i.e., the flux would extrapolate linearly to 0 at a distance d beyond the boundary. This zero flux boundary condition is more straightforward and is can be written mathematically as:

If d is not negligible, physical dimensions of the reactor are increased by d and extrapolated boundary is formulated with dimension R_e = R + d and this condition can be written as Φ(R + d) =  Φ(R_e) = 0.

It may seem the flux goes to 0 at an extrapolated length beyond the boundary. This interpretation is not correct. The flux cannot go to zero in a vacuum, because there are no absorbers to absorb the neutrons. The flux only appears to be heading to the zero value at the extrapolation point.

Note that, the equation d ≈ 0.7104 λ_tr is applicable to plane boundaries only. The formulas for curved boundaries can differ slightly, however, the difference is small unless the radius of curvature of the boundary is of the same order of magnitude as the extrapolated length.

Typical values of the extrapolated length:

The most common moderators have following diffusion coefficients (for thermal neutrons):

D(H₂O) = 0.142 cm

D(D₂O) = 0.84 cm

D(Be) = 0.416 cm

D(C) = 0.916 cm

The thermal neutron extrapolated lengths are given by:

d ≈ 0.7104 λ_tr = 0.7104 x 3 x D

therefore:

H₂O: d ≈ 0.30 cm

D₂O: d ≈ 1.79 cm

Be: d ≈ 0.88 cm

C: d ≈ 1.95 cm

As can be seen, this approximation is valid when the dimension L of the diffusing medium is much larger than the extrapolated length, L >> d.

This conditions are often used to eliminate unnecessary functions from solutions.

It must be added, as J must be continuous, the flux gradient will show a jump if the diffusion coefficients in both media differ from each other.

On this special interface we shall apply an albedo boundary condition to represent the neutron reflector. Albedo, the latin word for “whiteness”, was defined by Lambert as the fraction of the incident light reflected diffusely by a surface.

In reactor engineering, albedo, or the reflection coefficient, is defined as the ratio of exiting to entering neutrons and we can express it in terms of neutron currents as:

For sufficiently thick reflectors, it can be derived, that albedo becomes

where D_refl is the diffusion coefficient in the reflector and the L_refl is the diffusion length in the reflector.

If we are not interested in the neutron flux distribution in the reflector (let say in the slab B) but only in the effect of the reflector on the neutron flux distribution in the medium (let say in the slab A), the albedo of the reflector can be used as a boundary condition for the diffusion equation solution. This boundary condition is similar to the vacuum boundary condition, i.e. Φ(R_albedo) = 0, where R_albedo = R + d_e and

For x > 0, this diffusion equation has two possible solutions exp(x/L) and exp(-x/L), which give a general solution:

Φ(x) = Aexp(x/L) + Cexp(-x/L)

Note that, B is not usually used as a constant, because B is reserved for a buckling parameter. To determine the coefficients A and C we must apply boundary conditions.

From finite flux condition (0≤ Φ(x) < ∞), that required only reasonable values for the flux, it can be derived, that A must be equal to zero. The term exp(x/L) goes to ∞ as x ➝∞  and therefore cannot be part of a physically acceptable solution for x > 0. The physically acceptable solution for x > 0 must then be:

Φ(x) =  Ce^-x/L

where C is a constant that can be determined from source condition at x ➝0.

If S₀ is the source strength per unit area of the plane, then the number of neutrons crossing outwards per unit area in the positive x-direction must tend to S₀ /2 as x ➝0.

So that the solution may be written:

We can replace the 3D Laplacian by its one-dimensional spherical form, because there is no dependence on angle (whether polar or azimuthal). The source is assumed to be an isotropic source (there is the spherical symmetry).  The flux is then a function of radius – r only, and therefore the diffusion equation (outside the source) can be written as (everywhere except r = 0):

If we make the substitution Φ(r) = 1/r ψ(r), the equation simplifies to:

For r > 0, this differential equation has two possible solutions exp(r/L) and exp(-r/L), which give a general solution:

Note that, B is not usually used as a constant, because B is reserved for a buckling parameter. To determine the coefficients A and C we must apply boundary conditions.

To find constants A and C we use the identical procedure as in the case of infinite planar source. From finite flux condition (0≤ Φ(r) < ∞), that required only reasonable values for the flux, it can be derived, that A must be equal to zero. The term exp(r/L)/r goes to ∞ as r ➝∞  and therefore cannot be part of a physically acceptable solution for r > 0. The physically acceptable solution for r > 0 must then be:

Φ(r) =  Ce^-r/L/r

where C is a constant that can be determined from source condition at x ➝0.

If S₀ is the source strength, then the number of neutrons crossing a sphere outwards in the positive r-direction must tend to S₀ as r ➝0.

So that the solution may be written:

Since there is no dependence on angle Θ and z-coordinate, we can replace the 3D Laplacian by its one-dimensional form, and we can solve the problem only in radial direction. The source is assumed to be an isotropic source. Since the flux is a function of radius – r only, the diffusion equation (outside the source) can be written as (everywhere except r = 0):

This differential equation is called the Bessel’s equation and it is well known to mathematicians. In this case, the solutions to the Bessel’s equation are called the modified Bessel functions (or occasionally the hyperbolic Bessel functions) of the first and second kind, I_α(x) and K_α(x) respectively.

For r > 0, this differential equation has two possible solutions I₀(r/L) and K₀(r/L), the modified Bessel functions of order zero, which give a general solution:

To find constants A and C we use the identical procedure as in the case of infinite planar source. From finite flux condition (0≤ Φ(r) < ∞), that required only reasonable values for the flux, it can be derived, that A must be equal to zero. The term I₀(r/L) goes to ∞ as r ➝∞  and therefore cannot be part of a physically acceptable solution for r > 0. The physically acceptable solution for r > 0 must then be:

Φ(r) =  C.K₀(r/L)

where C is a constant that can be determined from source condition at r ➝0.

If S₀ is the source strength, then the number of neutrons crossing a cylinder outwards in the positive r-direction must tend to S₀ as r ➝0.

So that the solution may be written:

where a is the real width of zone 1 and b the outer dimension of the diffusion environment including the extrapolated distance. With problems involving two different diffusion media, the interface boundary conditions play crucial role and must be satisfied:

At interfaces between two different diffusion media (such as between the reactor core and the neutron reflector), on physical grounds the neutron flux and the normal component of the neutron current must be continuous. In other words, φ and J are not allowed to show a jump.

1., 2. Interface Conditions

It must be added, as J must be continuous, the flux gradient will show a jump if the diffusion coefficients in both media differ from each other. Since the solution of these two diffusion equations requires four boundary conditions, we have to use two boundary conditions more.

3. Finite Flux Condition

The solution must be finite in those regions where the equation is valid, except perhaps at artificial singular points of a source distribution.  This boundary condition can be written mathematically as:

4. Source Condition

The presence of the neutron source can be used as a boundary condition, because it is necessary that all neutrons flowing through bounding area of the source must come from the neutron source. This boundary condition depends on the source geometry and for planar sourve can be written mathematically as:

For x > 0, these diffusion equations have the following appropriate solutions:

Φ₁(x) = A₁exp(x/L₁) + C₁exp(-x/L₁)

and

Φ₂(x) = A₂exp(x/L₂) + C₂exp(-x/L₂)

where the four constants must be determined with use of the four boundary conditions. The typical neutron flux distribution in a simple two-region diffusion problem is shown at the picture below.

Since the neutron current is equal to zero (J = -D∇Ф, where Ф is constant), the diffusion equation in the infinite uniform multiplying system must be:

The only solution of this equation is a trivial solution, i.e., Ф = 0, unless Σ_a = νΣ_f. This equation (Σ_a = νΣ_f) is known as the criticality condition for an infinite reactor and expresses the perfect balance (critical state) between neutron absorption and neutron production. This balance must be continuously maintained in order to have steady state neutron flux.

Infinite Multiplication Factor

In this section the infinite multiplication factor, k_∞, will be defined from another point of view than in section – Nuclear Chain Reaction.

As can be seen we can rewrite the diffusion equation in following way and we can define a new factor – k_∞ = νΣ_f / Σ_a:

A non-trivial solution of this equation is guaranteed when k_∞ = νΣ_f / Σ_a = 1. On the other hand we have no information about the neutron flux in such critical reactor. In fact the neutron flux can have any value and the critical uniform infinite reactor can operate at any power level. It should be noted this theory can be used for a reactor at low power levels, hence  “zero power criticality”.

In power reactor core, the power level does not influence the criticality of a reactor unless thermal reactivity feedbacks act (operation of a power reactor without reactivity feedbacks is between 10E-8% – 1% of rated power).

See also: Reactor Criticality 

The quantity B_g² is called the geometrical buckling of the reactor and depends only on the geometry. This term is derived from the notion that the neutron flux distribution is somehow ‘‘buckled’’ in a homogeneous finite reactor. It can be derived the geometrical buckling is the negative relative curvature of the neutron flux (B_g² = ∇²Ф(x) / Ф(x)). In a small reactor the neutron flux have more concave downward or ‘‘buckled’’ curvature (higher B_g²) than in a large one. This is a very important parameter and it will be discussed in following sections.

For x > 0, this diffusion equation has two possible solutions sin(B_gx) and cos(B_gx), which give a general solution:

Φ(x) = A.sin(B_gx) + C.cos(B_g x)

From finite flux condition (0≤ Φ(x) < ∞), that required only reasonable values for the flux, it can be derived, that A must be equal to zero. The term sin(B_gx) goes to negative values as x goes to negative values and therefore it cannot be part of a physically acceptable solution. The physically acceptable solution must then be:

Φ(x) =  C.cos(B_g x)

where B_g can be determined from vacuum boundary condition.

The vacuum boundary condition requires the relative neutron flux near the boundary to have a slope of -1/d, i.e., the flux would extrapolate linearly to 0 at a distance d beyond the boundary. This zero flux boundary condition is more straightforward and is can be written mathematically as:

If d is not negligible, physical dimensions of the reactor are increased by d and extrapolated boundary is formulated with dimension a_e/2 = a/2 + d and this condition can be written as Φ(a/2 + d) =  Φ(a_e/2) = 0.

Therefore, the solution must be Φ(a_e/2) =  C.cos(B_g .a_e/2) = 0 and the values of geometrical buckling, B_g, are limited to B_g = ^nπ/_{a_e}, where n is any odd integer. The only one physically acceptable odd integer is n=1, because higher values of n would give cosine functions which would become negative for some values of x. The solution of the diffusion equation is:

It must be added the constant C cannot be obtained from this diffusion equation, because this constant gives the absolute value of neutron flux. In fact the neutron flux can have any value and the critical reactor can operate at any power level. It should be noted the cosine flux shape is only in theoretical case in a uniform homogeneous reactor at low power levels (at “zero power criticality”).

In power reactor core (at full power operation), the neutron flux can reach, for example, about 3.11 x 10¹³ neutrons.cm^-2.s^-1, but this values depends significantly on many parameters (type of fuel, fuel burnup, fuel enrichment, position in fuel pattern, etc.).

The power level does not influence the criticality (k_eff) of a power reactor unless thermal reactivity feedbacks act (operation of a power reactor without reactivity feedbacks is between 10E-8% – 1% of rated power).

Let assume a uniform reactor (multiplying system) in the shape of a sphere of physical radius R. The spherical reactor is situated in spherical geometry at the origin of coordinates. In order to solve the diffusion equation, we have to replace the Laplacian by its spherical form:

We can replace the 3D Laplacian by its one-dimensional spherical form, because there is no dependence on angle (whether polar or azimuthal). The source term is assumed to be isotropic (there is the spherical symmetry). The flux is then a function of radius – r only, and therefore the diffusion equation can be written as:

The solution of the diffusion equation is based on a substitution Φ(r) = 1/r ψ(r), that leads to equation for ψ(r):

For r > 0, this differential equation has two possible solutions sin(B_gr) and cos(B_gr), which give a general solution:

From finite flux condition (0≤ Φ(r) < ∞), that required only reasonable values for the flux, it can be derived, that C must be equal to zero. The term cos(B_gr)/r goes to ∞ as r ➝0  and therefore cannot be part of a physically acceptable solution. The physically acceptable solution must then be:

Φ(r) =  A sin(B_gr)/r

The vacuum boundary condition requires the relative neutron flux near the boundary to have a slope of -1/d, i.e., the flux would extrapolate linearly to 0 at a distance d beyond the boundary. This zero flux boundary condition is more straightforward and is can be written mathematically as:

If d is not negligible, physical dimensions of the reactor are increased by d and extrapolated boundary is formulated with dimension R_e = R + d and this condition can be written as Φ(R + d) =  Φ(R_e) = 0.

Therefore, the solution must be Φ(R_e) =  A sin(B_gR_e)/R_e = 0 and the values of geometrical buckling, B_g, are limited to B_g = ^nπ/R_e, where n is any odd integer. The only one physically acceptable odd integer is n=1, because higher values of n would give sine functions which would become negative for some values of x before returning to 0 at R_e. The solution of the diffusion equation is:

It must be added the constant A cannot be obtained from this diffusion equation, because this constant gives the absolute value of neutron flux. In fact the neutron flux can have any value and the critical reactor can operate at any power level. It should be noted this flux shape is only in theoretical case in a uniform homogeneous spherical reactor at low power levels (at “zero power criticality”).

In power reactor core, the neutron flux can reach, for example, about 3.11 x 10¹³ neutrons.cm^-2.s^-1, but this values depends significantly on many parameters (type of fuel, fuel burnup, fuel enrichment, position in fuel pattern, etc.).

The power level does not influence the criticality (k_eff) of a power reactor unless thermal reactivity feedbacks act (operation of a power reactor without reactivity feedbacks is between 10E-8% – 1% of rated power).

Since there is no dependence on angle Θ, we can replace the 3D Laplacian by its two-dimensional form, and we can solve the problem in radial and axial directions. Since the flux is a function of radius – r and height – z only (Φ(r,z)), the diffusion equation can be written as:

The solution of this diffusion equation is based on use of the separation-of-variables technique, therefore:

where R(r) and Z(z) are functions to be determined. Substituting this into the diffusion equation and dividing by R(r)Z(z), we obtain:

Because the first term depends only on r, and the second only on z, both terms must be constants for the equation to have a solution. Suppose we take the constants to be v² and к², sum of these constants must be equal to B_g² = v² + к². Now we can separate variables, and the neutron flux must satisfy the following differential equations:

Solution for radial direction

The differential equation for radial direction is called the Bessel’s equation and it is well known to mathematicians. In this case, the solutions to the Bessel’s equation are called the Bessel functions of the first and second kind, J_α(x) and Y_α(x) respectively.

For r > 0, this differential equation has two possible solutions J₀(vr) and Y₀(vr), the Bessel functions of order zero, which give a general solution:

From finite flux condition (0≤ Φ(r) < ∞), that required only reasonable values for the flux, it can be derived, that C must be equal to zero. The term Y₀(vr) goes to -∞ as r ➝0  and therefore cannot be part of a physically acceptable solution. The physically acceptable solution for must then be:

R(r) =  AJ₀(vr)

The vacuum boundary condition requires the relative neutron flux near the boundary to have a slope of -1/d, i.e., the flux would extrapolate linearly to 0 at a distance d beyond the boundary. This zero flux boundary condition is more straightforward and is can be written mathematically as:

If d is not negligible, physical dimensions of the reactor are increased by d and extrapolated boundary is formulated with dimension R_e = R + d and this condition can be written as Φ(R + d) =  Φ(R_e) = 0.

Therefore, the solution must be R(R_e) =  A J₀(vR_e) = 0. The function of J₀(r) has several zeroes, the first is at r₁ = 2.405, and the second at r₂ = 5.6. However, because the neutron flux cannot have regions of negative values, the only physically acceptable value for v is 2.405/R_e. The solution of the diffusion equation is:

Solution for axial direction

The solution for axial direction have been solved in previous sections (Infinite Slab Reactor) and therefore it has the same solution. The solution of in axial direction is:

Solution for radial and axial directions

The full solution for the neutron flux distribution in the finite cylindrical reactor is therefore:

where B_g² is the total geometrical buckling.

It must be added the constants A and C cannot be obtained from the diffusion equation, because these constant give the absolute value of neutron flux. In fact the neutron flux can have any value and the critical reactor can operate at any power level. It should be noted this flux shape is only in theoretical case in a uniform homogeneous cylindrical reactor at low power levels (at “zero power criticality”).

In power reactor core, the neutron flux can reach, for example, about 3.11 x 10¹³ neutrons.cm^-2.s^-1, but this values depends significantly on many parameters (type of fuel, fuel burnup, fuel enrichment, position in fuel pattern, etc.).

The power level does not influence the criticality (k_eff) of a power reactor unless thermal reactivity feedbacks act (operation of a power reactor without reactivity feedbacks is between 10E-8% – 1% of rated power).