Mathematical Methods The Sturm-Liouville Problem Neutron Diffusion in a Nuclear Reactor

Mathematical Methods (10/24.539)

IX. The Sturm-Liouville Problem

Example 9.2 -- Neutron Diffusion in a Nuclear Reactor

Problem Description:

Determine the 1-group neutron flux distribution in an infinitely long homogeneous bare critical cylindrical reactor. Also show that the resultant eigenfunctions for this problem are orthogonal and find the norm of these orthogonal eigenfunctions.

Problem Solution:

Background

This example is from basic Nuclear Reactor Theory. In a nuclear reactor, ²³⁵U fissions by neutron bombardment and recoverable energy to drive a steam power plant is released. Neutrons are also released in the fission process, and these are used to initiate additional fission reactions. For steady state operation, the parasitic neutron losses in the system and the neutron absorption in the uranium fuel are balanced exactly by the neutron production rate from fission. Since it is the neutron distribution that determines the various interaction rates and the energy production rate, it becomes important when designing and operating a nuclear reactor to determine the neutron population throughout the system. One must solve a second order differential balance equation to obtain the desired neutron density or flux distribution. The simplest form of the neutron balance equation is the so-called Neutron Diffusion Equation.

One of the steady state critical ideal reactor geometries that can be treated via analytical means is a long bare cylindrical core model, as illustrated in Fig. 9.3. All the adjectives used to describe the system are needed to reduce the general, very complicated, particle balance equation into a form that can be treated analytically. For example, the phrase, steady state critical, indicates that the neutron population is constant in time and that the reaction is self-sustaining (i.e. the neutrons given off in the fission process are balanced perfectly so that they cause the same number of fissions in each generation, which produces the same number of new neutrons, and keeps the total neutron population constant - on the average). The term bare reactor means that the fueled core region is surrounded by a vacuum, giving an outer boundary condition of zero flux (the neutron flux is a measure of the neutron population). The long homogeneous description implies that the axial height is large relative to the radius and that the material properties are constant throughout the system. These conditions suggest that the variation of the neutron density in the axial and azimuthal directions is negligible, leaving a functional dependence on only one variable, or , where is the symbol used to represent the neutron flux. Also, symmetry in the system suggests that the neutron population will be the largest in the center of the reactor, which implies that the flux gradient is zero at r = 0.

Fig. 9.3 Basic geometry for the cylindrical bare critical reactor in Example 9.2.

The Neutron Flux Distribution

With all the above conditions and simplifications, the Neutron Diffusion Equation for a one energy group (assumes all the neutrons have the same average energy) bare critical reactor model becomes

where the or Laplacian operator for 1-D cylindrical geometry is given by

and B² is a constant that is related to material properties of the homogenous system.

Therefore, the particle balance equation for the 1-g 1-D neutron flux distribution is

with boundary conditions as shown in Fig. 9.3, or

In mathematical terms, this system is a 2^nd order homogeneous ODE with homogeneous boundary conditions. From our recent discussions, we recognize this as an eigenvalue or Sturm-Liouville problem!

To avoid any confusion with the above notation, let’s first convert this system into a standard Sturm-Liouville problem, using the notation from previous sections. In particular, if we let , then, after multiplication by r, the balance equation becomes

This is clearly a Sturm-Liouville problem with

and boundary conditions

We should also recognize this as an ordinary Bessel equation with . This can be seen more easily by rewriting the balance equation as

Therefore, the general solution for this problem can be written as

Also, we can compute the gradient, , as

where the standard derivative formulas for the ordinary Bessel functions have been used, or

As usual, we now apply the boundary conditions to the general solution. At r = x = 0, we have

where we have used the facts that J₁(x) approaches zero as and Y₁(x) approaches negative infinity as . Thus, the only way to satisfy this condition is to let . Thus, after applying only the symmetry condition, the solution profile reduces to

Now, applying the boundary condition at x = R gives

Since A₁= 0 would give a trivial solution, we must let to satisfy this condition. However, J_n(x) is an oscillatory function and it has an infinite number of zeros, which we will denote as for -- that is, represents values of x for which .

Therefore, from the above discussion, we have

as the eigencondition for this problem. This leads to

as the allowed eigenvalues (really ) with .

With the eigenvalues known [these can be approximated from a plot of J₀(x) or, more accurately, from a root finding algorithm applied to J₀(x)], the eigenfunctions are simply

Note that, for the real physical system, the neutron flux must be real and non-negative, and the only eigenfunction that is positive over the full domain, , is related to the fundamental mode (i.e. first eigenfunction). Since the first zero of the J₀(x) Bessel function occurs at x = 2.4048, the real neutron flux distribution in the physical system is

where A₁ is a normalization that is usually determined from the overall power level of the reactor. This latter expression is the desired solution to the first part of the problem description for Example 9.2 (see above).

Orthogonality of the Eigenfunctions

In addition to finding the physical neutron flux profile, from a mathematical viewpoint, it would also be nice to demonstrate that the eigenfunction solutions are indeed orthogonal functions with respect to the weight function p(x) = x.. If we take the general orthogonality relationship given in eqn. (9.7) for the general Sturm-Liouville problem and apply it to this problem, we have

The first term on the RHS of this expression vanishes because, by definition, for , and the second term is zero for two reasons - because of the zero coefficient and because . Therefore, since the full RHS vanishes, the ordinary Bessel functions are indeed orthogonal with respect to p(x) = x, or

Finally, let’s compute the normalization required when . In this case we need to evaluate the following integral:

To accomplish this we first obtain from the literature the following integral relationship involving the square of the J₀Bessel function:

and, for k = 0, this becomes

Now to evaluate the norm of the eigenfunction, we need to put the desired integral into this form. Therefore, we let

and substitution into the expression for the normalization gives

where the last equality uses the fact that . Thus, the desired eigenfunction normalization is given by

__________________________

This example completes our current study of the general Sturm-Liouville Problem. You will find that these ideas, especially those related to orthogonality and generalized Fourier series, will be very useful in lots of applications. They are also especially relevant for several techniques for solving PDEs. Thus, we will revisit some of these concepts in the next section on PDEs (the next section deals with the analytical solution of PDEs using the Separation of Variables method).

10/24.539 Lecture Notes by Dr. J. R. White, UMass-Lowell (updated November 1998).

Return to Online Courses