# Solving the 2-Group Diffusion Equation Numerically

Document Sample

```					    Solving the 2-Group Diffusion Equation Numerically

B. Rouben
McMaster University
EP 4D03/6D03
Nuclear Reactor Analysis
2008 Sept-Dec

2008 September                                    1
Contents
   We derive the finite-difference version of the
2-group diffusion equation and a method to
solve it numerically.

2008 September                                   2
Considering Real Reactors
    Up to now we have considered very simple cases of
homogeneous reactors.
    Real life is more complicated, and real reactors have
very complex, certainly not homogeneous, cores!
    Therefore, analytic solutions to the flux shape in real
reactors are not applicable. Instead, we must be
prepared to solve the diffusion equation numerically.
    Also, the 1-group treatment is not sufficiently
accurate for real reactors. We must do diffusion
calculations in at least 2 energy groups.

2008 September                                         3
Diffusion Equation in 2 Energy Groups
   Let’s make a small simplifying concession:
   We will assume that there is no upscattering from the thermal
group to the fast group, only downscattering (moderation)
from the fast group to the thermal group, and
   We will assume that all fissions are thermal fissions, i.e., are
induced by neutrons in the thermal group.
   Also, we will work in 2 spatial dimensions instead of 3.
   [These assumptions are not crucial, and do not simplify
the method significantly, they just simplify the picture a
bit.]

2008 September                                                     4
Notation
   Let’s adopt the following notation:
   Subscript 1  Fast group
   Subscript 2 Thermal group
   r  position in 2 spatial dimensions (x, y)
   D1(r), D2(r)  Group diffusion coefficients
   a1(r), a2(r)  Group absorption cross sections
   12(r)  Downscattering (moderation) cross section
   f(r)  Thermal fission cross section
   keff  Reactor multiplication constant

2008 September                                    5
Diffusion Equation in 2 Energy Groups
   The 2-group diffusion equation is then:
 f 2 r 2 r 
   D1 r 1 r   a1 r   12 r 1 r                          0
keff
   D2 r 2 r   a 2 r 2 r   12 r 1 r                      0 (1)

   We have to get a solution which satisfies these equations at all
points in the reactor.
   Now we need to have a reactor model. See a model for a small
“2-dimensional” reactor on the next slide. The reactor is
modelled as square (rather than near-cylindrical) for simplicity.

2008 September                                                                              6
Reactor Model

2008 September                   7
Finite-Difference Method
   We will apply the finite-difference method to solve for the flux at
the midpoints of the cells defining the reactor model. Note that
the assumption is made that the cells defined have homogeneous
properties. That is, all cells of the same colour have the same
properties; of course, outer-core, inner-core, and reflector
properties are different!
   The method should also produce the value of keff, which is not
known a priori.
   The finite-difference method gets its name from the fact that it
approximates derivatives by finite-difference ratios.
   The finite-difference method can be mesh-centred or edge-
centred, depending on whether it solves for the flux at the centre
or at the edge of the meshes. We will use the mesh-centred
method.

2008 September                                                 8
Finite-Difference Method
   The finite-difference method starts by integrating Eq. (1) over the volume of a
cell, let’s label it with superscript C:
                                              C 22C r 
cellC 
    D1C 1C r   C1  12 1C r   f
a
C
dr  0
                                                 keff      

    D  r               r   121C r dr  0
C C       C C         C
2 2       a2 2                                           ( 2)
cell C

   The integrals in Eq. (2) are of 2 kinds: those which involve the divergence
operator, and those which don’t.
   The integrals which do not involve the divergence are easy to deal with,
taking advantage of the homogeneous properties. For example:

        C22C r dr  C2 
a               a            2C r dr     (3)
cell C                        cell C

2008 September                                                                         9
Finite-Difference Method (cont’d)
   At this point we make the approximation that the volume
integral of the flux is equal to the value at the centre of the cell
times the volume of the cell, i.e.:
        2C r dr  2C  dr  2CV C   ( 4)
cell C                   cell C

where 2C  2C ( r  centre of cell C ), and
V C  Volume of cell C  x C  y C (5)

   See a sketch on the next slide of cell C and its immediate
neighbours on the left, right, top, and bottom (L, R, T, B
respectively). Those are the neighbours we will need to
calculate the leakage from cell C.

2008 September                                                    10
Cell & Neighbours

2008 September                       11
Integrals in the Divergence Operator
   Let us now deal with the integrals in the divergence operator,
i.e., integrals

    D  r dr                             D  r dr
C    C                               C C
1   1               and              2 2
cell C                                      cell C

   These are the integrals which give the fast and thermal leakage
out of cell C. Let’s work with the thermal group as an example.
The int egral     D2 2C r dr is the volume int egral of the divergence
C
cell C

of the current J over cell C. By Green ' s Theorem, this volume
int egral can be equated to a surface int egral of the dot product of

ˆ
J with the outward normal n to the surface :
     C C
       

     
    D2 2 dr     J dr   J  n dS (6).
ˆ
C                    C            Surface of C

This form of the int egral clearly shows it ' s the leakage out of C.
2008 September                                                                12
Leakage Calculation in Finite Difference
   Let’s do the surface integral by first doing the integral
over one portion of the surface – say the interface
between C and L.
   In order to calculate the leakage in finite difference, we
need to use the flux at the midpoint of each cell, and
also the flux at the interface.
   The basic figure we need for the calculation is shown in
the next slide.

2008 September                                         13
Cells C and L and Interface

2008 September                                 14
Leakage Calculation
   To calculate the leakage, we first calculate the current (in the x
direction) on each side of the interface CL, approximating the
derivative of the flux by the first-order difference in fluxes
between the midpoint of each cell and the (unknown) flux at the
interface. We equate these two currents, since the current must
be continuous at the interface:
Jx   Jx 
CL     CL

2CL  2L          2C  2CL
 D         L
 D
C
( 7)
x   2         L
x      2         C

2                 2
   We use this equation to solve for the interface flux. We get
D2 2L xC  D2 2C x L
L            C
2CL                                        (8)
2008 September
D2 xC  D2 x L
L        C
15
Leakage Calculation (cont’d)
   Substituting this back into the right-hand side of Eq. (7) to
calculate J (we could equally have used the left-hand side):
C
2 D2 C D2L2L x C  D2 2C x L 
C
Jx   C   2 
                              

x         D2 x  D2 x
L    C    C     L


  L C 2 C L  C2L 
2 DCD2L
( 9)
D2 x  D2 x
2

   To calculate the leakage from cell C across interface CL, we must
take -Jx (since the outward normal to CL is towards -x) and
multiply it by the “area” of the interface. This “area” is yC: the
model is 2-dimensional, so we can take all cell “widths” in the 3rd
dimension (z) to be unity. Therefore
2 D2 D2 y C
2C  2L  (10)
C  L
Leakage out of C across CL  L C
D2 x  D2 x L
C

2008 September                                                   16
Leakage Calculation (cont’d)
   We could write Eq. (10) as
Leakage out of C across CL  A2 2C  2L 
CL
(11)
2 D2 D2L y C
C
where we wrote A    CL
 L C                     (12)
D2 x  D2 x L
2                 C

   An exactly similar procedure can be used to derive the leakage
across the other interfaces, so that we can write in the end
Total Leakage out of C 
A2 2C  2L  A2 2C  2R   A2 2C  2T   A2 2C  2B 
CL               CR                CT                CB

 A2 2C  A2 2L  A2 2R  A2 2  A2 2B 
C         CL       CR       CT T    CB
(13)
where we wrote A2  A2  A2  A2  A2
C    CL   CR   CT   CB
(14)

2008 September                                                         17
Leakage Calculation (cont’d)
   If we put together everything we have derived, we can write the
finite-difference form of Eq. (1) as:
C 2V C2C
A1C C C1  12 V C1C 
1    a
C                 f
 A1CL1L  A1CR1R  A1CT 1T  A1CB1B
keff
 12V C1C 
C
A2 2C  C2V C2C  A2 2L  A2 2R  A2 2  A2 2CB (15 )
C
a
CL       CR       CT T    CB

   Eq. (15) must be solved for every (cell) point (C) in the reactor
model. This means that there is a system of [2 x the number of
points in the model] coupled equations.
   For the small 2-d model in slide 7, it’s a system of hundreds of
coupled equations, and for a realistic reactor model it can be
thousands or tens of thousands of coupled equations.
   It is impossible to solve such large systems in one shot (i.e.,
“invert the matrix of the problem” in one shot). Thus, the
solution must be by iteration.
2008 September                                                                          18
Iterative Solution
   The way that the iterative solution works is as follows:
    Make an initial guess for the flux values (e.g., take the flux value
to be the same at every point C).
    Scan through the reactor points C one by one, solving the 2x2
system in Eq. (15) for the 2 flux values at C, using on the right-
hand sides the current (latest) values of the flux at the neighbour
points L, R, T, B, then proceeding to the next point C. This
represents one “iteration”.
    Repeat the scan of the reactor for another iteration (using as
latest guess the flux values obtained at the previous iteration),
then another, etc., until “convergence”.
    Declare that convergence has been achieved when the ratio of
the flux in the current iteration to that in the previous iteration,
for every point in the reactor, is sufficiently close to 1 (i.e., the
difference from 1 is smaller than a desired convergence
criterion, e.g., 10-4 or 10-5.)
2008 September                                                   19
Points at External Reactor Boundaries
   Before we can actually embark upon the iterative
solution, there is one thing that we must do: deal with the
boundary points, which have no “neighbours”.
   To deal with the boundaries, we make use of the
extrapolation distance, dextr = 0.7104 tr (into a vacuum).
   Consider a point C at the left edge of the reactor. It has
no neighbour L, but we can pretend that it has a
neighbour L where the flux is 0 at the extrapolation
distance. For this “dummy neighbour cell”, the width
will then be x = 2*dextr, so that the midpoint is at dextr.
   See the sketch on the next slide.

2008 September                                         20
Dummy Neighbour Cell Outside Reactor
   We can pretend that the
point C at the boundary has a
neighbour cell outside the
reactor, with diffusion
coefficient equal to the one
in cell C, cell width 2*dextr,
and flux = 0 at the midpoint.
   i.e., for this dummy cell L:
DiL  DiC , i  1, 2
x L  2 * 0.7104tr
 6 * 0.7104DiC
( will depend on i )
iL  0
2008 September                       21
Back to Iterative Solution
   We can now perform the iterations, if we make and use a ring of
“dummy” neighbours outside the external boundary, with
properties as shown in the previous slide.
   [Note: we do not solve for the flux in the dummy cells, we know
that it is 0, we use them only as neighbours, as needed.]
   One more thing: We need to determine keff as we iterate. We can
start with a guess (e.g., keff = 1), and after each iteration j update
keff by the ratio of total flux in the reactor to that in the previous
iteration, as the flux is expected to change according to the change
in keff :
  1C , j  2C , j 
j
keff  keff 1
j               all C

                                  (16)
1
C , j 1

C , j 1
2
all C

   Note: It is normally found that keff converges much faster than the
individual flux values.
2008 September                                                                22
Summary
   The equations derived in this learning object are the
basis for the numerical solution of the diffusion
equation in 2 energy groups by the finite-difference
technique.
   We did make it slightly simpler by not having an
explicit fast-fission cross section or an upscattering
cross section.
   Do you think that the addition of such cross sections
would make the method very different?

2008 September                                           23
Course Project
   Write a FORTRAN program to solve the diffusion equation for
the 2-d reactor model shown in slide 7.
   The nuclear properties to be used are shown in the following
slide.
   The course project involves designing and writing the computer
program, writing the project report which explains your
computer program and gives your results in terms of the
converged flux values at each point and the keff value. You
should also report on how many iterations the program took to
converge.
   Note: Even though the mesh spacings in the model are all equal
(except for the dummy cells around the boundary of the model),
write your program so that it allows for non-uniform mesh
spacings.

2008 September                                             24
Model Properties for Course Project
    Use the following values for the nuclear properties in the
various types of cells (materials) in the reactor model in the
course project:

Core Region             D1      D2       a1        a2       12        f
Reflector          1.314   0.876     0.0      0.793E-4    0.0101      0.0
Outer Core             1.269   0.937   0.766E-3   0.406E-2   0.740E-2   0.476E-2
Inner Core             1.269   0.937   0.766E-3   0.409E-2   0.742E-2   0.472E-2

2008 September                                                       25
END

2008 September         26

```
DOCUMENT INFO
Shared By:
Categories:
Tags:
Stats:
 views: 28 posted: 1/4/2012 language: English pages: 26
How are you planning on using Docstoc?