Ứng dụng ngôn ngữ lập trình hướng đối tượng python cho bài toán mô phỏng sự truyền neutron trong lò hạt nhân trong điều kiện tới hạn với cấu trúc hình học phẳng

J. Sci. & Devel. 2015, Vol. 13, No. 6: 1016-1027 
Tạp chí Khoa học và Phát triển 2015, tập 13, số 6: 1016-1027 
www.vnua.edu.vn 
1016 
APPLICATION OF PYTHON PROGRAMMING TOOLS FOR CRITICALITY SIMULATION 
OF NEUTRON TRANSPORT IN NUCLEAR REACTOR WITH SLAB GEOMETRY 
Luong Minh Quan, Nguyen Thi Thanh 
Faculty of Information Technology, Viet Nam National University of Agriculture 
Email: lmquan83@gmail.com/thanhnt@vnua.edu.vn 
Received date: 22.07.2015 Accepted date: 03.09.2015 
ABSTRACT 
Monte Carlo criticality calculations use the power iteration method to determine the eigenvalue (݇௘௙௙) and 
eigenfunction (fission source distribution) of the fundamental mode. However, the main problems of this method are 
the slow convergence of fission source distribution from the initial guess to the stationary solution, and the correlation 
between successive cycles which results in an under-prediction bias in the confident intervals of the estimated 
response. In this paper, we presented the Wielandt's method aiming to accelerate the convergence of the Monte 
Carlo power iteration. The object-oriented programming called Python prototype, was used to describe the standard 
Monte Carlo criticality power iterations for mono-kinetic particles and to compare the results obtained by the two 
different methods of acceleration mentioned above. The Wielandt's method successfully suppressed the auto-
correlation, even though no gain in the figure of merit seemed to occur. 
Keywords: Eigenvalue, eigenfunction, Monte Carlo criticality, power iteration, Wielandt’s. 
Ứng dụng ngôn ngữ lập trình hướng đối tượng python 
cho bài toán mô phỏng sự truyền neutron trong lò hạt nhân 
trong điều kiện tới hạn với cấu trúc hình học phẳng 
TÓM TẮT 
Sử dụng phương pháp lặp theo hàm mũ để giải phương trình truyền neutron trong lò phản ứng hạt nhân với 
cấu trúc hình học phẳng, một chiều bằng phương pháp mô phỏng Monte Carlo để xác định trị riêng và hàm riêng 
(tương ứng với sự phân bố của nguồn phân hạch) tồn tại hai vấn đề chính là: (1) sự hội tụ chậm của nguồn phân 
hạch dự đoán ban đầu về phân bố thực; (2) mối liên hệ giữa các chu kỳ trong mỗi mô phỏng sẽ dẫn đến sai số khi 
đánh giá độ tin cậy của kết quả thu được. Do đó, nghiên cứu này đã đề cập tới một phương pháp mới, Wielandt, với 
mục đích tăng tốc quá trình hội tụ của nguồn ban đầu và đưa ra phương án mới để đánh giá mối liên hệ giữa các 
chu kì liên tiếp trong một mô phỏng độc lập. Ngôn ngữ lập trình hướng đối tượng có tên gọi Python đã được sử dụng 
để mô tả sự truyền của hạt neutron đơn năng, đồng thời so sánh kết quả thu được bằng hai phương pháp khác 
nhau. Kết quả cho thấy, phương pháp mới đã giải quyết thành công các vấn đề đặt ra của bài toán Monte Carlo. Tuy 
nhiên, bài toán về hệ số phẩm chất vẫn cần tiếp tục nghiên cứu thêm. 
Từ khóa: Bài toán tới hạn Monte Carlo, hàm riêng, lặp theo hàm số mũ, phương pháp Wielandt, trị riêng. 
1. INTRODUCTION 
Basically, the behaviour of nuclear reactor 
can be simulated by coupling neutron transport 
theory and thermal-hydraulics. The neutron 
transport theory studies the interaction of 
neutron with matter, which is described by the 
Boltzmann equation (Reuss, 2008). Therefore, 
the computer codes developed for nuclear 
reactor simulations have to solve this equation 
in order to accurately describe the neutron 
transport in the reactor. Nowadays, the 
development of powerful computer has 
enhanced the capability of not only 
Luong Minh Quan, Nguyen Thi Thanh 
1017 
deterministic codes to numerically solve the 
neutron transport equation, but also of 
stochastic codes with very detailed geometry 
and continuous energy description of the 
neutron collisions. These stochastic codes are 
known as Monte Carlo codes. 
The basic principle of Monte Carlo method 
in particle transport is to simulate the entire 
life of each particle from its initial emission 
until its death either by absorption or leakage 
from the system boundaries. The frequency, 
nature and outcome of every interaction that 
may occur during the particle life are randomly 
sampled according to algorithms derived from 
particle physics laws. When the process is 
repeated for a large number of particles, the 
averages of the obtained results yield a detailed 
description of the transport process. However, 
the main problem in Monte Carlo criticality 
calculations is the convergence of the power 
iteration (https://docs.python.org/2/tutorial/) and 
the technique used to converge an arbitrary 
initial guess source distribution to thecritical 
source distribution. After the source is 
converged, it is possible to tally desired 
quantities like the dominant eigenvalue of the 
Boltzmann equation݇௘௙௙, the corresponding 
eigenvector (the stationary flux distribution), 
reaction rates, etc. 
This convergence is known to be 
problematic when the dominante ratio (ratio 
between the second and the first eigenvalues of 
the system) is close to one. This undesired 
phenomenon occurs for large systems of the size 
of about 100 times the neutron mean free path. 
The consequence of the power iteration is, on 
one hand, the correlation between successive 
cycles, due to the fact that the fission source 
bank of the previous cycle is used as a source for 
current cycle. This can severely affect the 
variance estimation, and always leads to 
underestimate of the variance. On the other 
hand, there is a transient phase of the 
simulation necessary to converge the source 
distribution to the stationary distribution; 
during this phase the quantities of interest 
cannot be tallied and this corresponds to 
“inactive” cycles. This transient will be 
evaluated by performing the criticality 
simulations with two different initial source 
guess: a uniform shape and a cosine shape, 
which almost corresponds to the critical shape 
for our numerical test. 
Another method, named the Wielandt’s 
method, has been proposed in order to improve 
the slow convergence of the power iteration 
implement. This acceleration method along with 
the Monte Carlo standard are analysed in terms 
of the advantages and disadvantages of each in 
simple but meaningful test cases. 
2. BASIC THEORY OF NEUTRON TRANSPORT 
The transport equation describes the 
statistical behaviour of a large population of 
particles. The exact number of particles per unit 
volume is continuously varying with time, even at 
steady state conditions. Under steady state 
conditions, the number density of particles 
oscillates around an average value corresponding 
to the solution of the steady state transport 
equation. A solution of the transport equation is 
required in many fields of nuclear engineering, 
notably in reactor physics, insafety and criticality, 
and in radiation shielding and protection. 
The fundamental assumptions in neutron 
transport calculation are as following: neutrons 
ca nbe treated as point-like particles travelling 
along straight lines between two collision 
points, and all neutron-neutron interactions can 
be ignored. These assumptions lead to a linear 
transport equation. Additional usual 
assumptions often made are medium 
homogeneity by regions and time-independent 
conditions. 
2.1. The steady-state Boltzmann neutron 
transport equation 
The stationary behavior of neutral particles 
transported through matter is described by the 
linear steady-state transport Boltzmann 
equation: 
Application of Python Programming Tools for Criticality Simulation of Neutron Transport in Nuclear Reactor 
with Slab Geometry 
1018 
ߗ෠ .∇Ψ൫r, E, Ω෡൯ + Σ୲(r, E)Ψ൫r, E, Ω෡൯ 	=
	∫ ∫ Σୱ(r, E′ ⇢ E, Ω. Ω′)Ψ ቀr, E, Ω෡ ′ቁ dE′∞଴ dଶΩ′ସπ +Q(r, E, Ω෡) (1) 
where: 
- Ψ൫r, E, Ω෡൯ is the angular flux in six-
dimensional phase space 
- ܮΨ	 = 	 ߗ෠ .∇Ψ൫r, E,Ω෡൯, the streaming term or 
leakage term (describe the neutrons crossing 
the boundaries) 
- ܶΨ	 = 	 Σ௧(ݎ,ܧ)Ψ൫r, E,Ω෡൯, the collision term 
- 
ܵΨ	 = 	 ∫ ∫ Σୱ(r, Eᇱ ⇢ E,Ω.Ωᇱ)Ψ൫r, E,Ω෡ᇱ൯dEᇱஶ଴ dଶΩᇱ,ସ஠ 
the scattering term (describe the change in 
energy and direction of neutron particles) 
- MΨ	 = 	Q൫r, E,Ω෡൯, the neutron source term 
- Σୱ : differential scattering macroscopic 
cross section 
- Σ୲ : total macroscopic cross section 
For criticality problem, neutron fission 
source comes from the fissile materials during 
fission events: Q൫r, E, Ω෡൯ 	=
ଵ
୩౛౜౜
χ(୉)
ସπ
∬ υΣ୤(r, E′) Ψ ቀr, E′, Ω෡ ′ቁ dE′dଶΩ′ (2) 
The quantities of interest are: 
- the scalar flux distribution or 
eigenfunction: ߶(ݔ,ܧ) 
- the largest eigenvalue or multiplication 
factor: ݇௘௙௙ 
- various response rates: ܴ(ݔ) 	=
	∫ Σ(x, E)Φ(x, E)dE 
2.2. The solution of theneutron transport 
equation in special case of critical slab 
geometry 
In the case of one-dimensional 
homogeneous critical slab within one-group 
diffusion theory, the general form of Eq.(1) can 
be simplified to the simple form: 
ௗమథ(௫)
ௗ௫మ
+ ߯ଶ߶(ݔ) 	= 	0 (3) 
where ߯ଶ 	= 	 ൫జΣ೑ିΣೌ൯
஽
	= 	 (௞∞ିଵ)
௅మ
 is the material 
buckling. In general, ݇∞ 	= 	ߝ݌݂ߟ - the infinitive 
multiplication factor, in which: 
- ߝ: fast fission factor; 
- ݌: resonance escape probability 
- ݂: thermal utilization factor; 
- ߟ	 = 	 ቀజΣ೑
Σ౗
ቁ
௙௨௘௟
: the production factor 
By taking into account the condition of the 
vanishing of the flux at the boundaries and the 
symmetry of the problem with respect to the 
plan at x = 0; also the critical condition: ߯	 = 	 గ
ଶ௔
, 
the solution of the Eq.(3) is written as following: 
߶(ݔ) 	= 	 గ
ଶ௅
cos	(గ௫
ଶ௅
) (4) 
Fig. 1. The distribution of flux inside the critical slab 
Luong Minh Quan, Nguyen Thi Thanh 
1019 
2.3. The Power iteration method 
The standard Monte Carlo criticality 
calculation uses the fission source iteration that 
is a simple power iteration to find the 
fundamental (biggest) eigenvalue k0 and its 
associated eigenvector Ψ଴. 
Re-writing (1) by using abbreviations: (ܮ + ܶ − ܵ)Ψ	 = 	 ଵ
୩౛౜౜
MΨ (5) 
The principle of the power iteration method 
is to evaluate the neutron flux at iteration 
n+1,Ψ(୬ାଵ) by applying the operator F = (L+T–
S)-1M on the flux at iteration n,Ψ(୬). By 
supposing that the values of flux and ݇௘௙௙ 	are 
known from iteration n, we can find out the 
values for iteration n+1 according to the 
following procedures: 
Ψ(୬ାଵ) 	= 	 ଵ
୩౛౜౜
(౤) FΨ(୬) 
and kୣ୤୤(୬ାଵ) 	= 	kୣ୤୤(୬) ∫ୢ୚	୑Ψ(౤శభ)∫ୢ୚	୑Ψ(౤) (6) 
where V is the six-dimensional phase space 
(space, energy and direction). 
Supposing {ui } is the set of eigenvectors of 
the operator F and that it forms a basis, any 
initial flux distribution can be expanded in 
function of{ui }: 
Ψ(଴) 	= 	 ∑ a୧u୧ with a୧ 	= 	 ∫ ܸ݀	Ψ(଴)u୧ (7) 
Successive applications of the operator F to 
the initial distribution will converge to the 
fundamental eigenvalue and eigenvector. It can 
be shown that both flux and multiplication 
factor can be approximated by using [4]: 
	Ψ(௡ାଵ) ≈ ܥଵ ൤ݑ଴ + ௔భ௔బ ቀ௞భ௞బቁ(௡ାଵ) ∗ ݑଵ൨ and 
݇௘௙௙
(௡ାଵ) ≈ ݇଴ ൤1 + ܥଶ ௔భ௔బ ቀ௞భ௞బቁ(௡) ∗ ቀ௞భ௞బ − 1ቁ൨ (8) 
It is well-known that the dominant error 
terms for k0 decay faster than the dominant 
error terms for Ψ(଴), especially when the 
dominance ratio ρ = k1/k0 is close to 1, that is for 
very large system with characteristic 
dimensions of several (~100) times the neutron 
mean free path. For this reason, the 
fundamental eigenvalue for a high dominance 
ratio problem will appear to converge faster 
than the associated eigenvector, which will be 
more contaminated by the higher harmonics. 
3. MONTE CARLO SIMULATIONS 
3.1. The random walk for a single neutron 
The random walk process provides a 
faithful simulation of the behaviour of a single 
particle, from source to death. It is a continuous 
line made of straight free paths connected at 
collision locations (scattering, absorption, 
fission) of the domain and by the cross-section 
data corresponding to materials present in the 
domain. During this walk, the particle state is 
characterized by its position, direction of travel 
and energy. 
The random walks can be generated by: 
- sampling the travelled distance in space 
- sampling the type of reaction 
- sampling the direction of the particle after 
the collision from the continuous distribution 
function related to the sampled reaction 
Fig. 2a. Three modes of interaction 
of neutron 
3.2. Criticality calculations 
Nuclear criticality is the ability to sustain a 
chain reaction with fission neutrons. In this 
case, the Monte Carlo simulation proceeds in 
cycles, in which the source in each cycle is given 
by the fission neutron distribution calculated in 
the previous cycle. Each cycle comprises the 
simulation of a batch of neutrons containing N 
source fission neutrons each of them executted 
by a single random walk. The number of 
neutrons generated at the end of each cycle is 
generally not equal to the number ofs ource 
neutrons which started at the beginning of the 
cycle. The multiplication factor in cycle n is 
defined as ratio: 
fission 
absorption 
Application of Python Programming Tools for Criticality Simulation of Neutron Transport in Nuclear Reactor 
with Slab Geometry 
1020 
Fig. 2b. Random walk for a single neutron from the source to death 
݇௘௙௙ 	= 	 ௡௨௠௕௘௥	௢௙	௙௜௦௦௜௢௡	௡௘௨௧௥௢௡	௜௡	௚௘௡௘௥௔௧௜௢௡	ேାଵ௡௨௠௕௘௥	௢௙	௙௜௦௦௜௢௡	௡௘௨௧௥௢௡	௜௡	௚௘௡௘௥௔௧௜௢௡	ே (9) 
The quality of the initial spatial 
distribution of source neutrons is an important 
issue. A very poor initial source guess can cause 
several first cycles estimate of k eff to be 
extremely poor. This situation can occur when 
only a small fraction of the fission source has a 
chance to make fissions. 
4. ACCELERATION METHODS 
4.1. Standard Monte Carlo method 
The standard Monte Carlo method for 
criticality calculation uses the power iteration 
algorithm to converge the user defined source to 
the stationary solution. The basis idea of this 
standard method is to follow neutron particles, 
from the source to the end of its life by leakage 
from outer boundaries or absorption, etc, 
generation by generation. This is achieved by 
defining an initial fission bank from the user 
source definition. All the neutrons of the bank 
are followed one by one until all the particles 
are exhausted, and whenever a fission occurs 
the corresponding fission neutrons are put in 
the next cycle fission bank. This is described in 
the Fig.3(a,b): 
- All fission neutrons of the same 
generation are put in a cycle 
- Fission neutrons of the current cycle 
become the fission source of the next cycle 
iteration 
- Eigenvalue and eigenvector will be 
estimated at the end of each cycle. 
All neutron 
from the fission source sites are
 exhausted
Calculate: eigenvalue
 eigenvector
 Store fission 
source sites for 
using in the
 next cycle
Go to the next cycle
Inherit the fission source sites from
the previous cycle
Get a fission source data from the
fission source bank and start a particle
Follow this particle
Determine if the 
particle is killed
yes
no
no
yes
no
yes
Start a new cycle
Determine if fission occurs
(13) not satisfy
Luong Minh Quan, Nguyen Thi Thanh 
10 ... or calculation in the 1950s and 1960s. 
More details of this method were described by 
Brian et al. (2008). The main principles of 
criticality simulation by Wielandt's method is 
described in Fig.4: 
cycle 1 cycle 2 cycle 3
initial
guess
Set properties of the neutron source:
energy: E = 1.
position (x, 0, 0)
direction (isotropic in LAB)
volume index
Compute the end position of free path
Reaction ?
Compute the free path using l
Compare l and d
l <d
yes
scatteringabsorption branching
Determine the type of reaction
Find the new direction 
of secondary particles
Deactivate particle
(n,2n)
l >d particle still
 inside Slab? void
no
yes
d: is the distance
 to the boundary
Determine:
 volume index
 medium
MC Source
݈ = −ln	(ߦ)
Σ௧
Application of Python Programming Tools for Criticality Simulation of Neutron Transport in Nuclear Reactor 
with Slab Geometry 
1022 
Fig. 4. Compare MC and Wielandt’s of implementation 
The difference between standard Monte 
Carlo method and Wielandt's method lies in the 
way we treat the neutrons produced by fission. 
For MC, each fission neutron is put in the next 
cycle, but for the Wielandt's method, a fraction 
of these neutrons are kept in the current cycle 
and continued to be tracked in a single random 
walk by the probability condition: 
௖ܰ௨௥௥௘௡௧ 	= 	 ఔΣ೑Σ೟ ൬ ଵ௞೐೑೑ − ଵ௞೐൰ + ߦ ≤ 0,8 (10) 
݇௘ 	= 	 ݇௘௙௙ + 4ߪ(ߪ(݇௘௙௙)); ߪ is the variance of 
variance of keff 
If fission is determined to be occurring, the 
number of fission neutrons Ncurrent are tracked by 
the same random walks in the current cycle just 
like other source particles. These neutrons 
continue producing their progenies, and 
continue tracking within the current cycle until 
their death by leakage, absorption, etc, (Fig. 4). 
The chain of current-generation neutrons 
followed by such away will be treated as a part 
of the same initial history in the iteration, so 
that correlation effects between successive 
fissions will decrease. This improved statistical 
treatment should reduce the under-prediction 
bias in confidence intervals for criticality 
calculations. The value on the right hand side of 
Eq.(10) depends on the size of the slab, the 
larger the slab is the higher the value is. 
5. RESULTS AND DISCUSSION 
5.1. Numerical parameters 
5.1.1. Properties of homogeneous slab 
To simplify our problem, we consider a 1D 
homogeneous slab with void boundary conditions. 
In our criticality simulations, the particles are 
supposed to be mono-kinetic (MKParticle: its 
velocity is always constant and it doesn't change 
with scattering collisions), and the fission and 
scattering interactions are supposed isotropic. The 
homogeneous medium was chosen with the 
following properties in Tab. 1. 
5.1.2. Parameters of simulations 
- Simulations will be performed for three 
different lengths of slab: L = 10; 50; 100 (cm) 
- Two different shapes of initial fission 
source guess (uniform and cosine shape) are 
also investigated to analyze the behavior of the 
transient phase. 
Tab. 1. Macroscopic cross section for homogeneous slab 
Types of reaction total absorption scattering branching ̅ߥ 
Macroscopic cross section (cm-1) Σ௧ 	= 	1,0 Σ௔ 	= 	0,2 Σ௦ 	= 	0,8 Σ௙ 	= 	0,1 2 
Standard Monte Carlo method
Start from a fission source site
Killed
(2)
(1)
(0)
Wielandt's method
Start from a fission source site
Killed(0)
(1) (0)
(0)
(1)
(1)
(0)
Killed
Killed
Luong Minh Quan, Nguyen Thi Thanh 
1023 
5.2. Numerical results 
5.2.1. The dependence of keff on the size of 
the critical slab 
By definition of multiplication factor, this 
quantity is the ratio of the total number of 
particles produced by fission at the end of the 
cycle to the number of initial particles; this ratio 
is strongly sensitive to the effect of leakage from 
the boundaries. For the same slab material 
properties, the smaller the length of the slab, 
the bigger the effect of leakage, and finally the 
lower the value of keff. We also note that the 
smaller the size of the slab (in mean free paths), 
the better the particle can explore the entire 
slab in a single generation. For large slabs, it 
takes several generations for a particle to travel 
from one side of the system to the other. Note 
that: 1pcm = 10-5 
The adoption of the Wielandt's method 
leads to an increase of statistics because of the 
tracking of additional particles in each cycle. 
Consequently, the fluctuations and the 
standard deviation of multiplication factor are 
less than the one implemented by conventional 
Monte Carlo simulations. Those results could be 
seen clearly on Fig. 5 for the case L = 100 cm. 
5.2.2. Figure of Merit (FOM) 
In Monte Carlo simulations, the variance 
associated with the estimated quantities is 
(asymptotically) linearly proportional to the 
inverse of the number of histories 
simulated:	ߪଶ 	= 	1
√ܰ
ൗ . On the other hand, the 
simulation time is roughly proportional to the 
number of histories. This has led, for the 
comparison of different simulations and/or 
different methods in terms of efficiency, to the 
definition of a criterion called Figure of Merit 
(FOM):	ܨܱܯ	 = 	 ଵ
ఙమ்
 (11), in which T is the 
simulation time. 
When comparing two set of simulations on 
the same problem, the run with a bigger FOM is 
considered more efficient (smaller variance 
and/or smaller computing time). To evaluate the 
efficiency of our different implemented methods, 
we have carried out a group of 1000 
independent simulations, 1000 particles, 1000 
cycles for slab length L = 50 cm. 
The FOM for MC and Wielandt's are not 
much different and the gains in variance due to 
the improved statistics are compensated by the 
extra time spent in the simulations (Table 3). 
Tab. 2. The values of keff for different methods 
Slab length (L in cm) MC (ߪ	݅݊	݌ܿ݉) Wielandt (ߪ	݅݊	݌ܿ݉) 
10 0,88931 ± 3,2 0,88845 ± 1,5 
50 0,99371 ± 20,1 0,99298 ± 5,6 
100 0,99832 ± 28,9 0,99779 ± 16,6 
Fig. 5. The fluctuation keff in MC and Wielandt for L = 100 cm 
Application of Python Programming Tools 
with Slab Geometry 
1024 
Tab. 3. FOM of 1000 independent simulationsfor L = 50 cm
Methods
Time (h)
ߪଶ(10-10
FOM 
5.2.3. Auto-correlation estimation
The cycle to cycle correlation is one of the 
main problems of the power iteration method, 
due to the fact that the source of neutrons for 
the next cycle comes from neutron which have 
suffered fission in the current cycle. This 
correlation depends on the size of the system: 
more cycles are required for neutron to be 
propagated from one side to the other of larger 
systems. This correlation length (
2012) can be estimated at the end of the Monte 
Carlo simulation by introducing an empirical
correlation function as following: 
ܥ௞ 	= 	 〈(Φ೔ିఓ)∗(Φ೔శೖିఓ)〉ఙమ (12)
where: 
-	Φ௜ the flux at cycle i; 
- ߤ	the expected value of variable 
-	ߪଶ the variance of variable Φ
- ݇	 ≥ 1 
A low value of the auto
indicates that the sequence made by consecutive 
values of	Φ௜ 	with i = 1,2 ,... , k
Fig. 6. Correlation vs k for bin 11 (start 
curve) and bin 41 (diamon curve) f
for Criticality Simulation of Neutron Transport in Nuclear Reactor
 MC Wielandt’s 
 100 419 
) 404 87 
69 77 
Malvagi et al., 
Φ 
; 
-correlation 
 are weakly 
correlated. We see that the auto
decreases with the value of the lag k (Fig.6). 
When this value is below the cutoff threshold of 
0.2, we empirically consider that is 
uncorrelated. On Fig.7 we plot the value of the 
lag for which the auto-correlation descends 
under the cutoff versus the bin number 
corresponding to the slab with the length of 100 
and 50 cm, respectively, simulated by standard 
MC and Wielandt's method. This figure is 
consistent with the theory discussed above: the 
larger the slab, the higher the value of 
correlation length. By adopting Wielandt's 
method to treat the cycle
this quantity strongly decreases from 86 
for L = 100 cm (high DR problem) and from 40 
to 10 for L = 50 cm. From our simulations for L 
= 10 cm, there is no correlation problem for any 
of the methods. It also remarks that the 
correlations are at the lowest near the 
boundaries of the system (t
effects of the boundary conditions) and at its 
center (this is where the second harmonics has 
a node). 
or L = 100 
Fig. 7. Correlation cut
as function of number of bin cells
-correlation 
-to-cyclecorrelation, 
to 45 
his is due to the 
-off (0,2) 
Luong Minh Quan, Nguyen Thi Thanh 
1025 
5.2.4. Wielandt's method implementation to 
cancel the auto-correlation 
The basic idea of Wielandt's method is to 
keep neutrons to stay longer in each cycle and, 
thus, spreading out more in space compared to 
the standard MC. To be better understood the 
effects of this method on the correlation lengths, 
we carried out the simulations for L = 100 cm 
(large size) for a set of 100 independent 
simulations, varying the fraction of fission 
neutrons kept in each cycle. 
Fig. 8 plotts the max of correlation cutoff in 
the slab as a function of the fraction of fission 
neutrons kept in each cycle. We can conclude 
that for high dominant ratio systems (larger 
size), the Wielandt's method is a useful way to 
decrease the cycle by cycle correlations. 
5.3. Fission source convergence 
5.3.1. Characteristics of the fission source 
convergence as a function of the initial 
source guess and slab size 
The power iteration technique allows to 
converge a user defined source to the critical 
source distribution in Monte Carlo simulation. 
For 1D homogeneous slab geometry, void 
boundary condition and for both uniform and 
cosine initial source guess distribution, we 
expect the equilibrium flux distribution to 
assume a cosine form in space (predicted by 
diffusion theory). The convergence of the 
criticality simulation is estimated by measuring 
the flux in each spatial bin. In this part, the 
behaviour of fission source convergence 
depending on the slab length and initial source 
guess is investigated. The results are described 
in Fig. 9. 
The Fig. 9(a) shows the flux convergence and 
fluctuations, for L = 50 cm, for the both uniform 
and cosine shape implemented by standard MC 
and Wielandt's for simulations. These indicate 
the importance of choosing the right initial 
source guess for MC to quickly converge the flux 
distribution. However, by adopting Wielandt’s 
method, we can find the robust initial guess of 
user defined source because of increasing the 
statistics of each cycle iteration leading to 
increase the spatial distribution. 
For the case of L = 100 cm, in Fig.9(b), the 
results are presented for the same uniform 
source guess but for the set of 1 and 100 
independent simulations. The fluctuations and, 
specially, transient phase are hightly enhanced. 
We obtain very large size system with a high 
dominance ratio. 
Fig. 8. Correlation length for cut-off factor 0,2 and L = 100 cm 
0% 10% 20% 30% 40% 50% 60% 70% 80% 90% 100%
0
10
20
30
40
50
60
70
80
90
100
co
rr
el
at
io
n 
le
ng
th
 fo
r L
=1
00
 c
m
Application of Python Programming Tools for Criticality Simulation of Neutron Transport in Nuclear Reactor 
with Slab Geometry 
1026 
L = 50 cm Monte Carlo Wielandt’s 
Uniform 
(1 simu) 
(a) 
Cosine 
(1 simu) 
(a) 
L = 100 
cm 
uniform 
(1 simu) 
(b) 
uniform 
(100 
simus) 
(b) 
Fig. 9. The distribution of flux in special grid as function of cycle iteration 
For each independent simulation, standard 
MC results in the source distribution are very 
concentrated in some small region of the slab. 
This behavior can be explained as follow: in the 
regions where more fission neutrons exist, there 
will be a high probability for the neutrons 
interacting with medium to produce new fission 
neutrons. On the contrary, by adopting 
Wielandt's method for the same initial number 
of particles: moreparticles will be tracked in one 
cycles; these particles will expand in a greater 
spatial region of the slab system. Consequently, 
the fission source distribution will spread out 
more in a given number of cycles. 
6. CONCLUSION AND FUTURE OUTLOOK 
The main objective of this article was to use 
a Python prototype to implement the standard 
Monte Carlo criticality power iterations for 
mono-kinetic particles and to compare with the 
Wielandt's method of accelerating the power 
iteration. Numerical simulations of a critical 
slab of varying lengths where then used to 
compare the performances of these methods. 
The author was successful to implement 
and modify the available source code in the 
laboratory of SERMA in CEA for the standard 
Monte Carlo which is the most basic method of 
Luong Minh Quan, Nguyen Thi Thanh 
1027 
implementation to describe the interaction of 
neutrons in the core. However, due to the slow 
convergence of the source, especially in the case 
of the linear initial shape and large size of the 
geometry, it needs to accelerate the process of 
the convergence to get the better results in the 
same level of the competence. 
The Wieland’s method implemented was only 
concentrated in the direction of accelerating the 
source convergence and there were no ideas to 
estimate the multiplication factor, keff, which is 
one of the most important quantity to know if the 
new method is suitable or not. 
The three keys objectives has been obtained 
as following: 
- Based on the basic of Wielandt’s method 
to accelerate and to reduce the time of 
convergence of any kind of initial guess source. 
- Proposing the new way to estimate the 
value of keff inside the code and compare it to 
the one of MC 
- Estimating the correlation between cycle 
to cycle by increasing the statistics for the same 
initial condition with the old method of MC. 
There is a need of future work for a better 
idea of the biases of those methods that can 
introduce into the simulation as a function of their 
respective parameters: fraction of fission neutrons 
kept in the current cycle for Wielandt's. Another 
question not explored is the optimisation between 
the number of cycles run in a simulation and the 
number of independent simulations in order to 
minimize the final variance, which will certainly 
depend on the method. 
REFERENCES 
Brian, C. Kiedrowski, Forrest B (2008). Brown. Using 
Wielandt’s Method to Eliminate Confidence 
Interval Underprediction Bias in MCNP5 
Criticality Calculations. 
Brown F. (2005). Fundamentals of Monte Carlo 
Particle Transport. LAUR -05-4983, Los Alamos 
National Laboratory. 
Fausto Malvagi, Eric Dumonteil, Francois-Xavier 
Hugot (2012). Les bonnes pratiques dans les 
calculs critiques en Monte Carlo. Commissariat a 
l’Energie Atomique DEN/DANS/DM2S/SERMA. 
https://docs.python.org/2/tutorial/ 
https://root.cern.ch/drupal/ 
Kitada T., T. Takeda (2000). Effective Convergence of 
Fission Source Distribution in Monte Carlo 
Simulation. J. Nucl. Sci. Techlo., 38(5): 324 - 329. 
Paul Reuss (2008). Neutron Physics. Institute National 
des Sciences et Techniques Nuclear.