## Solutions of the Diffusion Equation – Multiplying Systems

In previous section it has been considered that the environment is **non-multiplying**. In non-multiplying environment neutrons are emitted by a neutron source situated in the center of coordinate system and then they freely diffuse through media. We are now prepared to consider **neutron diffusion** in **multiplying system**, which contains fissionable nuclei (i.e. **Σ**_{f }**≠ 0**). In this multiplying system we will also study spatial distribution of neutrons, but in contrast to non-multiplying environment these neutrons can trigger **nuclear fission reaction**.

**stable, self-sustained fission chain reaction**in a multiplying system (in a nuclear reactor) is that

**exactly every fission initiate another fission**. The minimum condition is for each nucleus undergoing fission to produce, on the average, at least one neutron that causes fission of another nucleus. Also the number of fissions occurring per unit time (the reaction rate) within the system must be constant.

This condition can be expressed conveniently in terms of **the multiplication factor**. The infinite multiplication factor is the ratio of the **neutrons produced by fission** in one neutron generation to the number of **neutrons lost through absorption** in the preceding neutron generation. This can be expressed mathematically as shown below.

It is obvious** the infinite multiplication factor** in a multiplying system is a measure of the change in the fission neutron population from one neutron generation to the subsequent generation.

**k**. If the multiplication factor for a multiplying system is_{∞}< 1**less than 1.0**, then the**number of neutrons is decreasing**in time (with the mean generation time) and the chain reaction will never be self-sustaining. This condition is known as**the subcritical state**.

**k**. If the multiplication factor for a multiplying system is_{∞}= 1**equal to 1.0**, then there is**no change in neutron population**in time and the chain reaction will be**self-sustaining**. This condition is known as**the critical state**.

**k**. If the multiplication factor for a multiplying system is_{∞}> 1**greater than 1.0**, then the multiplying system produces**more neutrons**than are needed to be self-sustaining. The number of neutrons is exponentially increasing in time (with the mean generation time). This condition is known as**the supercritical state**.

In this section, we will solve the following **diffusion equation**

in various geometries that satisfy the **boundary conditions**. In this equation** ν** is number of neutrons emitted in fission and **Σ _{f}** is macroscopic cross-section of fission reaction.

**Ф**denotes a

**reaction rate**. For example a fission of

^{235}U by thermal neutron yields

**2.43 neutrons.**

It must be noted that we will solve the diffusion equation without any external source. This is very important, because such equation is a **linear homogeneous equation** in the flux. Therefore if we find one solution of the equation, then any multiple is also a solution. This means that the** absolute value** of the neutron flux **cannot possibly be deduced** from the diffusion equation. This is totally different from problems with external sources, which determine the absolute value of the neutron flux.

We will start with simple systems (planar) and increase complexity gradually. The most important assumption is that all neutrons are lumped into a **single energy group**, they are emitted and diffuse at** thermal energy** (**0.025 eV**). Solutions of diffusion equations in this case provides an illustrative insights, how can be the neutron flux distributed in a reactor core.

It is known the fission neutrons are of importance in any chain-reacting system. Neutrons trigger the nuclear fission of some nuclei (^{235}U, ^{238}U or even ^{232}Th). What is crucial the fission of such nuclei produces **2, 3 or more** free neutrons.

But not all neutrons are released **at the same time following fission**. Even the nature of creation of these neutrons is different. From this point of view we usually divide the fission neutrons into two following groups:

**Prompt Neutrons.**Prompt neutrons are emitted**directly from fission**and they are emitted within**very short time of about 10**.^{-14}second**Delayed Neutrons.**Delayed neutrons are emitted by**neutron rich fission fragments**that are called**the delayed neutron precursors**. These precursors usually undergo beta decay but a small fraction of them are excited enough to undergo**neutron emission**. The fact the neutron is produced via this type of decay and this happens**orders of magnitude later**compared to the emission of the prompt neutrons, plays an extremely important role in the control of the reactor.

**Radiative Capture.**Most absorption reactions result in the loss of a neutron coupled with the production of one or more gamma rays. This is referred to as a**capture reaction**, and it is denoted by**σ**._{γ}

**Neutron-induced Fission Reaction.**Some nuclei (fissionable nuclei) may undergo a fission event, leading to two or more fission fragments (nuclei of intermediate atomic weight) and a few neutrons. In a fissionable material, the neutron may simply be captured, or it may cause nuclear fission. For fissionable materials we thus divide the absorption cross section as**σ**._{a}= σ_{γ}+ σ_{f}

**uniform infinite reactor**, i.e. uniform infinite multiplying system without an external neutron source. This system is in a Cartesian coordinate system and under these assumptions (no neutron leakage, no changes in diffusion parameters) the

**neutron flux**must be

**inherently constant**throughout space.

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

**uniform reactor**(multiplying system) in the shape of a slab of physical

**width**

**a**in the x-direction and infinite in the y- and z-directions. This reactor is situated in the center at x=0. In this geometry the flux does not vary in y and z allowing us to eliminate the y and z derivatives from ∇

^{2}. The flux is then a

**function of x only**, and therefore the Laplacian and diffusion equation can be written as:

The quantity **B _{g}^{2}** 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**). In a small reactor the

_{g}^{2}= ∇^{2}Ф(x) / Ф(x)**neutron flux**have more concave downward or ‘‘buckled’’ curvature (

**higher B**) than in a large one. This is a very important parameter and it will be discussed in following sections.

_{g}^{2}For x > 0, this diffusion equation has two possible solutions **sin(B _{g}x)** and

**cos(B**, which give a general solution:

_{g}x)**Φ(x) = A.sin(B _{g}x) + 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 _{g}x)** 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) = 0Therefore, 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**, where n is any

_{g}=^{nπ}/_{a_e}**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**^{13 }**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 _{g}r)** and

**cos(B**, which give a general solution:

_{g}r)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 _{g}r)/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 _{g}r)/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}) = 0Therefore, the solution must be **Φ(R _{e}) = A sin(B_{g}R_{e})/R_{e} = 0** and the values of geometrical buckling, B

_{g}, are limited to

**B**, where n is any

_{g}=^{nπ}/R_{e}**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**^{13 }**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).

**uniform reactor**(multiplying system) in the shape of a cylinder of physical radius

**R and height H.**This finite cylindrical reactor is situated in cylindrical geometry at the origin of coordinates. In order to solve the

**diffusion equation**, we have to replace the Laplacian by its cylindrical form:

**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 ^{2}** and

**к**, sum of these constants must be equal to

^{2}**B**. Now we can

_{g}^{2}= v^{2}+ к^{2}**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 _{0}(vr)** and

**Y**, the Bessel functions of order zero, which give a general solution:

_{0}(vr)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 _{0}(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 _{0}(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}) = 0Therefore, the solution must be **R(R _{e}) = A J_{0}(vR_{e}) = 0**. The function of

**J**has several zeroes, the first is at

_{0}(r)**r**, and the second at r

_{1}= 2.405_{2}= 5.6. However, because the neutron flux cannot have regions of negative values, the only physically acceptable value for

**v**is

**2.405/R**. The solution of the diffusion equation is:

_{e}**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}^{2} **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**^{13 }**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).

**Nuclear and Reactor Physics:**

- J. R. Lamarsh, Introduction to Nuclear Reactor Theory, 2nd ed., Addison-Wesley, Reading, MA (1983).
- J. R. Lamarsh, A. J. Baratta, Introduction to Nuclear Engineering, 3d ed., Prentice-Hall, 2001, ISBN: 0-201-82498-1.
- W. M. Stacey, Nuclear Reactor Physics, John Wiley & Sons, 2001, ISBN: 0- 471-39127-1.
- Glasstone, Sesonske. Nuclear Reactor Engineering: Reactor Systems Engineering, Springer; 4th edition, 1994, ISBN: 978-0412985317
- W.S.C. Williams. Nuclear and Particle Physics. Clarendon Press; 1 edition, 1991, ISBN: 978-0198520467
- G.R.Keepin. Physics of Nuclear Kinetics. Addison-Wesley Pub. Co; 1st edition, 1965
- Robert Reed Burn, Introduction to Nuclear Reactor Operation, 1988.
- U.S. Department of Energy, Nuclear Physics and Reactor Theory. DOE Fundamentals Handbook, Volume 1 and 2. January 1993.

**Advanced Reactor Physics:**

- K. O. Ott, W. A. Bezella, Introductory Nuclear Reactor Statics, American Nuclear Society, Revised edition (1989), 1989, ISBN: 0-894-48033-2.
- K. O. Ott, R. J. Neuhold, Introductory Nuclear Reactor Dynamics, American Nuclear Society, 1985, ISBN: 0-894-48029-4.
- D. L. Hetrick, Dynamics of Nuclear Reactors, American Nuclear Society, 1993, ISBN: 0-894-48453-2.
- E. E. Lewis, W. F. Miller, Computational Methods of Neutron Transport, American Nuclear Society, 1993, ISBN: 0-894-48452-4.