NRA Assignment.

Department of Nuclear Engineering.
Assignment Group # 17.

“Numerical Solution for Neutron Diffusion
Equation for a One-dimensional Non-
homogenous reactor.”
Subject: Nuclear Reactor Analysis.
Submitted by:
Fazal e Rabbi (MS- Process Engineering).
Hamza Imtiaz Butt (MS- Process Engineering).
8-23-2022
1
Derivation of Difference Equation.
INTRODUCTION:
The one-speed neutron diffusion equation for a non-homogeneous reactor in three dimensions
can be written as:
1 𝜕𝜙
− 𝛻𝐷𝛻𝜙 + 𝛴𝑎𝜙 = 𝑆(𝑥, 𝑦, 𝑧) -------------------------- (A)
𝑣 𝜕𝑡
Due to the non-homogeneous nature of the reactor, the diffusion co-efficient (D) and the absorption
cross-section (Σa) are functions of space and hence the second term cannot be easily resolved. First
of all, we need to simplify our case by taking two assumptions:
1. Steady state i.e the time derivative of flux = 0.
2. One-dimensional case i.e the spatial dependence terms can be taken for x-direction only.
Applying these assumptions, equation (A) can be reduced to:
𝑑 𝑑
− 𝐷(𝑥). 𝜙 (𝑥) + 𝛴𝑎(𝑥). 𝜙(𝑥) = 𝑆(𝑥) --------------------- (B)
𝑑𝑥 𝑑𝑥
In order to determine the neutron flux as a function of x, we will solve the above differential
equation using numerical techniques.
DISCRETIZATION:
We begin by dividing the one-dimensional space into N+1 number of mesh points ranging from
xo to xN. The mesh intervals between consecutive points refer to individual mesh sizes which are
variable as we are considering the non-homogeneous case where the fuel, moderator and absorber
have different mean free paths and hence the mesh size will vary from point to point. Larger the
mean free path, larger the mesh size can be selected as there will be lesser probability of
interactions inside fine mesh. This discretization is shown in the following diagram:
∆1 ∆2 ∆i ∆N
xo x1 x2 xi-1 xi xN-1 xN
2
Whereas;
xi = Mesh point where the values of S, D, Σ are known or will have to be calculated.
∆i = Mesh size or interval between two consecutive mesh points.
We have to take into account the fact that for fuel and absorber materials, fine mesh size is required
to capture variation in flux. Hence, we further divide the meshes and define:
∆i/2 ∆i/2 ∆i+1/2 ∆i+1/2
xi-1 xi - ∆i/2 xi xi + ∆i+1/2 xi+1
SOLUTION:
Now we integrate equation (B) for the required mesh i.e from xi - ∆i/2 to xi + ∆i+1/2.
xi + ∆i+1/2
𝑑 𝑑
∫ [− {𝐷(𝑥). 𝜙 (𝑥)} + 𝛴𝑎(𝑥). 𝜙(𝑥) = 𝑆(𝑥)] 𝑑𝑥 (C)
𝑑𝑥 𝑑𝑥
xi − ∆i/2
The right-hand side of the equation:

xi + ∆i+1/2 xi + ∆i+1/2
∫ 𝑆(𝑥)𝑑𝑥 = 𝑆(𝑥𝑖) ∫ 𝑑𝑥
xi − ∆i/2 xi − ∆i/2
As previously mentioned, the values of the parameters S, D, Σ and ϕ are known for the mesh
points i, i-1 and i+1. Hence, S(xi) = Si and the solution of the above integral is as follows:
xi + ∆i+1/2
𝑥𝑖 + ∆(𝑖 + 1) − 𝑥𝑖 + ∆𝑖
∫ 𝑆(𝑥)𝑑𝑥 = 𝑆𝑖 [ ]
2
xi − ∆i/2
3
xi + ∆i+1/2
∆(𝑖 + 1) + ∆𝑖
∫ 𝑆(𝑥)𝑑𝑥 = 𝑆𝑖 [ ]
2
xi − ∆i/2
Similarly, the second term on the left-hand side becomes:

xi + ∆i+1/2
∆(𝑖 + 1) + ∆𝑖
∫ 𝛴𝑎(𝑥). 𝜙(𝑥)𝑑𝑥 = 𝛴𝑎𝑖 𝜙𝑖 [ ]
2
xi − ∆i/2
The first term is a little complex as it contains two derivatives. The negative sign has to be
neglected here. It can be resolved as follows:
xi + ∆i+1/2
𝑑 𝑑𝜙 𝑑𝜙 xi + ∆i+1/2
∫ (−𝑑𝑥) {𝐷(𝑥) } = |𝐷(𝑥) | (D)
𝑑𝑥 𝑑𝑥 𝑑𝑥 xi − ∆i/2
xi − ∆i/2
In equation (D), the derivative term can be solved as:
𝑑𝜙 xi + ∆i+1/2 𝑑𝜙 𝑑𝜙
|𝐷(𝑥) | = |𝐷(𝑥) | − |𝐷(𝑥) |
𝑑𝑥 xi − ∆i/2 𝑑𝑥 xi + ∆i+1/2 𝑑𝑥 xi− ∆i/2
Now based on the first principle and the definition of the derivative, for small interval (mesh size),
we can solve the two derivatives as:
𝑑𝜙 𝜙(𝑖 + 1) − 𝜙𝑖
| | =
𝑑𝑥 xi + ∆i+1/2 ∆(𝑖 + 1)
𝑑𝜙 𝜙𝑖 − 𝜙(𝑖 − 1)
| | =
𝑑𝑥 xi− ∆i/2 ∆𝑖
Where ϕi, ϕi-1 and ϕi+1 are known values of flux at xi, xi-1 and xi+1 respectively. Similarly for known
values for diffusion coefficient (Di, Di-1 and Di+1) at these points, we can assume an average
midpoint value for D(x) and write:
𝐷𝑖 + 𝐷(𝑖 + 1)
|𝐷(𝑥)|xi + ∆i+1/2 = = 𝐷𝑖, 𝑖 + 1
2
𝐷𝑖 + 𝐷(𝑖 − 1)
|𝐷(𝑥)|xi− ∆i/2 = = 𝐷𝑖, 𝑖 − 1
2
4
Hence, in equation (D), the integral becomes:

xi + ∆i+1/2
𝑑 𝑑𝜙 𝜙(𝑖 + 1) − 𝜙𝑖 𝜙𝑖 − 𝜙(𝑖 − 1)
∫ (−𝑑𝑥) {𝐷(𝑥) } = (𝐷𝑖, 𝑖 + 1) [ ] − (𝐷𝑖, 𝑖 − 1) [ ]
𝑑𝑥 𝑑𝑥 ∆ ( 𝑖 + 1) ∆𝑖
xi − ∆i/2
This equation can then be expanded as:

xi + ∆i+1/2
𝑑 𝑑𝜙 (𝐷𝑖, 𝑖 − 1) 𝐷𝑖, 𝑖 + 1 𝐷𝑖, 𝑖 − 1
∫ (−𝑑𝑥) {𝐷(𝑥) } = 𝜙(𝑖 − 1) − ( + ) 𝜙𝑖
𝑑𝑥 𝑑𝑥 ∆𝑖 ∆(𝑖 + 1) ∆𝑖
xi − ∆i/2
(𝐷𝑖, 𝑖 + 1)
+ 𝜙(𝑖 + 1) (E)
∆(𝑖 + 1)
Combining all of these results into equation (C), we can write:

𝑎(𝑖, 𝑖 − 1)𝜙(𝑖 − 1) + 𝑎(𝑖, 𝑖)𝜙𝑖 + 𝑎(𝑖, 𝑖 + 1)𝜙(𝑖 + 1) (F)
Where;
𝐷𝑖 + 𝐷(𝑖 − 1) 1
𝑎(𝑖, 𝑖 − 1) = − ( )
∆𝑖 ∆𝑖 + ∆(𝑖 + 1)
𝐷(𝑖 + 1) + 𝐷𝑖 𝐷𝑖 + 𝐷(𝑖 − 1) 1
𝑎(𝑖, 𝑖) = 𝛴𝑖 + ( + )
∆(𝑖 + 1) ∆𝑖 ∆𝑖 + ∆(𝑖 + 1)
𝐷(𝑖 + 1) + 𝐷𝑖 1
𝑎(𝑖, 𝑖 + 1) = − ( )
∆𝑖 ∆𝑖 + ∆(𝑖 + 1)
Now we have obtained a set of ‘N-1’ difference equations for a system of ‘N+1’ fluxes. Hence,
we need two more equations to solve this system. For this purpose, we can use the vacuum
extrapolated boundary condition while carefully assuming that the mesh points xO (ϕO = 0) and xN
(ϕN = 0) are placed at the extrapolated boundary. A more general boundary condition can be
considered such as non-reentrant current, for which the final two equations become:
𝑎0,0 𝜙0 + 𝑎0,1 𝜙1 = 𝑆0 and
𝑎𝑁,𝑁−1 𝜙𝑁−1 + 𝑎𝑁,𝑁 𝜙𝑁 = 𝑆𝑁 .

5
RESULTS AND DISCUSSION:

Hence, with the number of equations equal to the number of variables, we can consider our system
closed. Such sets of difference equations appear often as a result of one-dimensional steady state
diffusion equations i.e second order differential equations. The coefficient ‘a’ depends mainly on
the scheme used to derive them. If the mesh interval is small, these differences will be negligible
in actual calculations. Since the spatial variation of flux is generally characterized by diffusion
length (L). We can choose a mesh interval less than L for practical purposes. Similar difference
equations can also be derived for curvilinear geometries using cylindrical or spherical coordinates.

NRA Assignment.

Uploaded by

Copyright:

Available Formats

NRA Assignment.

Uploaded by

Document Information

Original Description:

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

NRA Assignment.

Uploaded by

Copyright:

Available Formats

Department of Nuclear Engineering.

Assignment Group # 17.

∆i/2 ∆i/2 ∆i+1/2 ∆i+1/2

xi-1 xi - ∆i/2 xi xi + ∆i+1/2 xi+1

The right-hand side of the equation:

Similarly, the second term on the left-hand side becomes:

In equation (D), the derivative term can be solved as:

Hence, in equation (D), the integral becomes:

This equation can then be expanded as:

Combining all of these results into equation (C), we can write:

𝑎𝑁,𝑁−1 𝜙𝑁−1 + 𝑎𝑁,𝑁 𝜙𝑁 = 𝑆𝑁 .

RESULTS AND DISCUSSION:

You might also like