COURSE EVALUATION

• Your chance to tell us what you think about the course

• Information about the evaluation process

• Link to the evaluation platform

• Deadline: October 13th

1
PROJECT

• 6 distinct exercises need to be carried out by each student to produce a portfolio of scripts (and associated write-up < 5 pages each)
solving standard reactor physics problems
1. Attenuation of neutrons from a 1D planar source
2. 1D/1G reactor slab Week Date
5 08-Oct
3. 1D/2G reactor slab 7 30-Oct
4. 1D/2G reactor slab with a loading pattern 10 12-Nov
10 20-Nov
5. Fuel Evolution with exposure in a homogeneous media 12 03-Dec
6. Reactor evolution during a cycle (space and time) 13 10-Dec

• Collaboration between students is highly encouraged but we expect an original script per student, e.g. no copy paste. Questions about
the coding exercises will be asked during the oral exam.

• For each exercise, the scripts (w/o write-up) are submitted through Moodle before the next one and a selection of them will be
checked by the course instructors

• LateX Templates are provided for each exercise

2
PROJECT

• The final portfolio of final scripts and associated write-ups should be delivered by the exam date (January 2025)
➢ 50% of the final grade
➢ Evaluation (equal weight for each exercise) based on
✓ The instructor can run the scripts independently 3 /6
✓ The script produces accurate results 1 /6
✓ The quality of the write-up (description of the script - input/algorithm/output - answers to questions) 2 /6

3
PROBLEM #1
Attenuation of a 1D planar source

• Consider a plate of diffusive material of finite thickness (between x=0 and x=x0). Neutrons, arriving from a uniform planar source,
located at x=0 and emitting S #.cm2s-1, diffusing within the plate.
• Σ𝑎 = 0.02 𝑐𝑚−1
• Σ𝑠 = 4 𝑐𝑚−1
• 𝑥0 = 10 𝑐𝑚 (physical boundary)
• 𝑆 = 1000 𝑛. 𝑐𝑚−2. 𝑠 −1 (4π emission)

1
• Diffusion coefficient D = 3 Σ𝑠+Σ𝑎
1
• Diffusion length L = 3 Σ𝑠 +Σ𝑎 Σ𝑎
2
• Extrapolated length d = 2D =
3 Σ𝑠 +Σ𝑎

4
PROBLEM #1
Learning outcome

1. Discretization of the diffusion equation with the finite difference method including expression of B.C.
2. Implementation of a numerical solver; verification against an analytical solution
3. Awareness of issues associated with mesh refinement (mainly in terms of accuracy)

5
FINITE DIFFERENCE METHOD

𝜕 𝜕
Diffusion Equation in 1D 𝐷 𝑥 𝛷 𝑥 − Σ𝑎 𝑥 𝛷 𝑥 + 𝑄 𝑥 = 0
𝜕𝑥 𝜕𝑥

Finite difference approximation

𝛷𝑖−1 𝛷𝑖 𝛷𝑖+1 𝜕 𝑓 𝑥 + 𝛥𝑥 − 𝑓 𝑥
𝑓 𝑥 ≈
𝜕𝑥 𝛥𝑥
𝛷𝑖−1ൗ 𝑥𝑖+1ൗ x
2
2

𝛥𝑥
FINITE DIFFERENCE METHOD
Streaming term at (i)

𝜕 𝜕
We want to determine 𝐷 𝑥 𝛷 𝑥
𝜕𝑥 𝜕𝑥
𝛷𝑖−1 𝛷𝑖 𝛷𝑖+1
Derivative at (i-1/2)
𝜕 𝛷𝑖 − 𝛷𝑖−1/2 𝛷𝑖−1ൗ 𝑥𝑖+1ൗ x
2
𝛷 𝑥 ቤ ≈ 2
𝜕𝑥 𝑖−1ൗ
𝛥𝑥/2
2
Derivative at (i+1/2)
𝜕 𝛷𝑖+1/2 − 𝛷𝑖 𝛥𝑥
𝛷 𝑥 ቤ ≈
𝜕𝑥 𝑖+1ൗ
𝛥𝑥/2
2

Streaming term at (i)

𝜕 𝜕 1 𝛷𝑖+1/2 − 𝛷𝑖 𝛷𝑖 − 𝛷𝑖−1/2
𝐷 𝑥 𝛷 𝑥 ≈ 𝐷𝑖 − 𝐷𝑖
𝜕𝑥 𝜕𝑥 𝛥𝑥 𝛥𝑥/2 𝛥𝑥/2

7
FINITE DIFFERENCE METHOD
Streaming term at (i)

Determination of 𝛷𝑖−1/2 through current continuity at i-1/2

𝛷𝑖 − 𝛷𝑖−1/2 𝛷𝑖−1/2 − 𝛷𝑖−1

−2𝐷𝑖 = −2𝐷𝑖−1
𝛥𝑥 𝛥𝑥

𝐷𝑖−1 𝛷𝑖−1 + 𝐷𝑖 𝛷𝑖
𝛷𝑖−1/2 =
𝐷𝑖−1 + 𝐷𝑖 𝛷𝑖−1 𝛷𝑖 𝛷𝑖+1
Similarly for 𝛷𝑖+1/2 𝛷𝑖−1ൗ 𝑥𝑖+1ൗ x
2
2

𝐷𝑖+1 𝛷𝑖+1 + 𝐷𝑖 𝛷𝑖
𝛷𝑖+1/2 =
𝐷𝑖+1 + 𝐷𝑖 β𝑖+1 𝛥𝑥
So it comes
𝐷𝑖+1 𝛷𝑖+1 + 𝐷𝑖 𝛷𝑖
𝛷𝑖+1/2 − 𝛷𝑖 − 𝛷𝑖 1 2𝐷 𝐷
𝐷𝑖+1 + 𝐷𝑖 𝑖 𝑖+1
𝐷𝑖 = 𝐷𝑖 = 𝛷 − 𝛷𝑖
𝛥𝑥/2 𝛥𝑥/2 𝛥𝑥 𝐷𝑖+1 + 𝐷𝑖 𝑖+1
β𝑖−1 β are surface properties, there
are 2 values per mesh point
𝛷𝑖 − 𝛷 1
𝑖−2 1 2𝐷𝑖 𝐷𝑖−1 β𝑖−1
−𝐷𝑖 =− 𝛷 − 𝛷𝑖−1 = − 𝛷𝑖 − 𝛷𝑖−1
𝛥𝑥/2 𝛥𝑥 𝐷𝑖−1 + 𝐷𝑖 𝑖 𝛥𝑥
8
FINITE DIFFERENCE METHOD
Discretized Form

Finally, the streaming term is: 𝜕 𝜕 1 𝛷𝑖+1/2 − 𝛷𝑖 𝛷𝑖 − 𝛷𝑖−1/2

𝐷 𝑥 𝛷 𝑥 ≈ 𝐷 − 𝐷𝑖
𝜕𝑥 𝜕𝑥 𝛥𝑥 𝑖 𝛥𝑥/2 𝛥𝑥/2
𝜕 𝜕 1
𝐷 𝑥 𝛷 𝑥 ≈ 2 β𝑖+1 𝛷𝑖+1 − 𝛷𝑖,𝑗 β𝑖+1 + β𝑖−1 + β𝑖−1 𝛷𝑖−1
𝜕𝑥 𝜕𝑥 𝛥𝑥

Back to the diffusion equation:

𝜕 𝜕
𝐷 𝑥 𝛷 𝑥 − Σ𝑎 𝑥 𝛷 𝑥 + 𝑄 𝑥 = 0
𝜕𝑥 𝜕𝑥
𝛷𝑖−1 𝛷𝑖 𝛷𝑖+1
So it comes
𝛷𝑖−1ൗ 𝑥𝑖+1ൗ x
2
2
1
β 𝛷 − 𝛷𝑖,𝑗 β𝑖+1 + β𝑖−1 + β𝑖−1 𝛷𝑖−1 − Σ𝑎,𝑖 𝛷𝑖 + 𝑄𝑖 = 0
𝛥𝑥 2 𝑖+1 𝑖+1

𝛥𝑥

Re-arranging the terms:

β𝑖+1 + β𝑖−1 β𝑖+1 β𝑖−1

𝛷𝑖 + Σ𝑎,𝑖 − 𝛷 − 𝛷 = 𝑄𝑖,𝑗
𝛥𝑥 2 𝛥𝑥 2 𝑖+1 𝛥𝑥 2 𝑖−1
9
FINITE DIFFERENCE METHOD
Boundary Conditions

Zero incoming partial current (vacuum) on the RHS

𝛷𝑖+1/2 𝛷𝑖+1/2 − 𝛷𝑖 1 𝐷𝑖
𝐽𝑥−𝑖+1ൗ 𝛷𝑖+1/2 = 𝛷
=
4
+ 𝐷𝑖
𝛥𝑥
=0 1 𝐷𝑖 𝛥𝑥 𝑖
2
4 + 𝛥𝑥

Fixed net current on the LHS 𝛷𝑖−1 𝛷𝑖 𝛷𝑖+1

𝛥𝑥 𝛷𝑖−1ൗ 𝑥𝑖+1ൗ x
𝛷𝑖 − 𝛷𝑖−1/2 2
𝐽𝑥𝑖−1ൗ = 𝐶 = −𝐷𝑖 𝛷𝑖−1/2 = 𝐶 + 𝛷𝑖 2
2 𝛥𝑥/2 2𝐷𝑖

𝛥𝑥

10
FINITE DIFFERENCE METHOD
B.C. at the RHS
Mesh formulation assuming zero incoming partial current at i+1/2
1𝐷𝑖 1
𝛷𝑖+1/2 = 𝛷 = 𝛷𝑖
1 𝐷𝑖 𝛥𝑥 𝑖 𝛥𝑥 1
4 + 𝛥𝑥 4𝐷𝑖 + 1 𝐵𝐶𝑖 =
𝛥𝑥
2 4𝐷 + 1
𝑖
Streaming term
𝛷𝑖+1/2 − 𝛷𝑖 𝐷𝑖,𝑗 1 1
𝐷𝑖 = 𝛷 −1 =− 𝛷𝑖
𝛥𝑥/2 𝛥𝑥/2 𝑖 𝛥𝑥 𝛥𝑥
4𝐷𝑖 + 1 2 4𝐷 + 1
𝑖 𝛷𝑖−1 𝛷𝑖 𝛷𝑖+1
𝛷𝑖 − 𝛷 1
𝑖−2 β𝑖−1
−𝐷𝑖 =− 𝛷𝑖 − 𝛷𝑖−1 𝛷𝑖−1ൗ 𝑥𝑖+1ൗ x
𝛥𝑥/2 𝛥𝑥 2
2

𝜕 𝜕 1 β𝑖−1
𝐷 𝑥 𝛷 𝑥 ≈ −𝐵𝐶𝑖 𝛷𝑖 − 𝛷𝑖 − 𝛷𝑖−1
𝜕𝑥 𝜕𝑥 𝛥𝑥 𝛥𝑥
𝛥𝑥
So it comes

𝐵𝐶𝑖 𝛥𝑥 + β𝑖−1 β𝑖−1

𝛷𝑖 + Σ𝑎,𝑖 − 𝛷 = 𝑄𝑖
𝛥𝑥 2 𝛥𝑥 2 𝑖−1
11
FINITE DIFFERENCE METHOD
B.C. at the LHS
Mesh formulation assuming fixed net current at i-1/2,j
𝛥𝑥/2
𝛷𝑖−1/2 = 𝐶 + 𝛷𝑖
𝐷𝑖

Streaming term

𝛷𝑖+1/2 − 𝛷𝑖 β𝑖+1
𝐷𝑖 = 𝛷𝑖+1 − 𝛷𝑖
𝛥𝑥/2 𝛥𝑥
𝛷𝑖−1 𝛷𝑖 𝛷𝑖+1
𝛷𝑖 − 𝛷 1
𝑖−2 𝑥𝑖+1ൗ x
−𝐷𝑖 =𝐶 𝛷𝑖−1ൗ 2
𝛥𝑥/2 2

𝜕 𝜕 1
𝐷 𝑥 𝛷 𝑥 ≈ 𝛥𝑥 2 β𝑖+1 𝛷𝑖+1 − 𝛷𝑖 + 𝐶𝛥𝑥
𝜕𝑥 𝜕𝑥
𝛥𝑥
So it comes

β𝑖+1 β𝑖+1 𝐶
𝛷𝑖 + Σ𝑎,𝑖 − 𝛷 = + 𝑄𝑖
𝛥𝑥 2 𝛥𝑥 2 𝑖+1 𝛥𝑥
12
SOLVING THE DISCRETIZED NEUTRON DIFFUSION EQUATION

a111 + a122 = S1
a211 + a222 + a233 = S2
a322 + a333 + a344 = S3
For everything there is to
 know about iterative methods:
check this
a N −1, N − 2 N − 2 + a N −1, N −1 N −1 = S N −1

In matrix form: 𝐀Φ = 𝐒 Φ = 𝐀−1 𝐒 Φ S

 a11 a12 Direct Inversion

 1   S1 
    
 2   S 2 
- Considerable computer power
 a21 a22 a23 - Computer round-off errors
    

a32 a33 a34 
 3  =  3  S
 a43 a44 a45 • JacobiIterative solutions:  4   S 4 
 • Gauss-Seidel
   
    
 A a N −1, N − 2 a N −1, N −1   N −1   S N −1 

PROBLEM #1
Attenuation of a 1D planar source

• Consider a plate of diffusive material of finite thickness (between x=0 and x=x0). Neutrons, arriving from a uniform planar source,
located at x=0 and emitting S #.cm2s-1, diffusing within the plate.
• Σ𝑎 = 0.02 𝑐𝑚−1
• Σ𝑠 = 4 𝑐𝑚−1
• 𝑥0 = 10 𝑐𝑚 (physical boundary)
• 𝑆 = 1000 𝑛. 𝑐𝑚−2. 𝑠 −1 (4π emission)

1
• Diffusion coefficient D = 3 Σ𝑠+Σ𝑎
1
• Diffusion length L = 3 Σ𝑠 +Σ𝑎 Σ𝑎
2
• Extrapolated length d = 2D =
3 Σ𝑠 +Σ𝑎

14
PROBLEM #1
Questions

• Find the analytical solution – plot it

➢ Report the value of the flux at x0 (scientific format with 4 significant digits)
➢ Report the value of the flux at 0 (scientific format with 4 significant digits)
• Discretize the diffusion equation using the finite difference method
• Express the relationship between 𝛷𝑖 , 𝛷𝑖+1 , 𝛷𝑖−1
➢ Compute the associated coefficients of the matrix A (provide 4 significant digits) for a mesh size of 0.1cm
• Express the boundary conditions at the source and at the physical boundary of the material
➢ Relationship between 𝛷𝑖 , 𝛷𝑖−1/2 and 𝛷𝑖 , 𝛷𝑖+1/2 respectively
➢ Compute the associated coefficients of the matrix A (provide 4 significant digits) for 𝛷𝑖 , 𝛷𝑖+1 , 𝛷𝑖−1 and a mesh size
of 0.1cm (100 mesh points)
• Write the associated numerical solver in the language of your choice; report its algorithm
➢ Use an existing solver to determine the solution of the linear system of equations

• Compare the analytical and numerical solutions

➢ Plot the difference between the solutions for a mesh size of 0.1 cm.
➢ Report the value of the flux at x0 (4 significant digits)
➢ Report the value of the flux at 0 (scientific format with 4 significant digits)
➢ Evolution with mesh size of the absolute error at 𝛷 𝑥0 . What do you notice ? How do you explain it?

15
HINTS

A_{i,j}: 1.6604e+01
A_{i-1,j}: -8.2919e+00
A_{i+1,j}: -8.2919e+00

ϕ𝑎𝑛𝑎 1.05 = 7.3299 103 ;

ϕ𝑛𝑢𝑚 1.05 = 7.3303 103 ;
𝑒(1.05) = 3.7901 10−1

16