Tutorial of general-purpose reactor physics code system CBZ
Go CHIBA
June 6, 2024

1

Introduction

When neutron collides with nucleus, nuclear ﬁssion reaction might occur with given probability and this nucleus isdivided into (generally) two fragments. These fragments, known as ﬁssion fragments, possess a huge amount of kinetic
energy of about 167 MeV, and nuclear power utilizes this energy ﬁnally as electric power. 1
Nuclear fuel in which ﬁssion reactions occur is generally contained in a pressure vessel. In a pressure vessel, variousequipments are located, and a small part around nuclear fuels is called a reactor core, and a reactor core is a “heart”
of a nuclear power plant.
In nuclear power plants, the number of neutrons which cause ﬁssion reactions is controlled to maintain ﬁssion
reactions. That is, to properly control ﬁssion reactions, we have to know and understand behavior of neutrons in areactor core. Neutron behavior in a reactor core is described by the Boltzman-type neutron transport equation, so wehave to solve this equation somehow. To solve this equation to know neutron behavior is called reactor core analysis orneutronics analysis, and these are one of important subjects in the ﬁeld of nuclear engineering.
Many code systems for reactor core analyses have been developed so far in the world. One of the most famouscode systems might be CASMO/SIMULATE developed by STUDVIK. This code system has been widely used for
actual light water power reactors analyses in many countries. In Japan, various code systems have been developed, andthe SRAC code of Japan Atomic Energy Agency (JAEA) and AEGIS/SCOPE-2 developed by Nuclear Fuel Industry,
Nuclear Engineering Ltd. and Nagoya university would representative.
The author worked at JAEA in the past. There he had developed a code system CBG by himself, and had publishedseveral papers with it. Based on the development experience of CBG, he began to develop a new code system CBZafter moving to Hokkaido University. This is the reason why both of the names “CBZ” and “CBG” are used throughthis tutorial. With CBZ, one can perform various kinds of reactor physics calculations such as criticality calculations,kinetic calculations, radiation shielding calculations and nuclear fuel burnup calculations. This tutorial is for users tolearn how to use the code system CBZ.

1 Strictly speaking, ﬁssion energy includes energy of γ-ray instantaneously emitted in the ﬁssion reaction and energy of β-ray and γ-ray
which are emitted from unstable ﬁssion fragments nuclides.

12Install

When one uses CBZ, he can do his calculations by using various modules provided in CBZ. This can be done for users
by writing a computer program with C++, and compiling and running it. In this compiling, a library ﬁle libCBG.a,which can be generated by compiling varous modules of CBZ, should be included. This is shown in Fig. 1. In thissection, we will compile the CBZ source program. How to compile computer programs written by users will be explainedin the following sections.

Fig. 1: Relation between the CBZ source ﬁles and computer programs prepared by users
After uncompressing and expanding the package named as cbz.xxxxyyzz.tgz, 2 you can get a directory CBZ whichcontains the following directories: 3
 CBG
 CBGCAL CBGLIB CBGLIB BURNIn the CBG directory, you can ﬁnd a src directory. In this directory, source codes of CBZ are located.
The CBGCAL directory is used to locate data to perform actual calculations, and in the distribution package you canﬁnd several directories such as Takeda.tutorial, lanl.tutorial and mza.tutorial.
The CBGLIB directory is used to store nuclear reaction data (nuclear data), CBZLIB 4 , and the CBGLIB BURN directory
is to store data ﬁles required in nuclear fuel burnup calculations.
Now let us compile CBZ. Please move to the CBG/src directory. 5 There is possibility that compilation has been
already done in the distribution package, so please type make clean to delete all the ﬁles previously generated. Then tostart compiling, please type make. After that, it starts compiling and a ﬁle libCBG.a might be generated after several
minutes. By typing ls *.a, you can conﬁrm that this ﬁle is successfully generated. When the CBZ source programs
are revised, users have to compile them again and newly generate libCBG.a.
This is the end of the install. Note that CBGLIB might not be included in the package because of its large size. Insuch a case, please contact with Chiba or others to get the CBGLIB data.

2 Please type tar

xvfz cbz.xxxxyyzz.tgz.

3 In the CBZ directory, many ﬁles are located in layer-structured directories.

When one wants to move all the data contained in the CBZ
directory, it is quite convenient if one can create one ﬁle including all these data. Doing “tar” is to summarize one directory including thedata into one ﬁle, and a summarized ﬁle is called a tar ﬁle. For example, when ones create a tar ﬁle for the CBZ directory, they have totype tar cvf cbz.tar CBZ, and then one ﬁle cbz.tar which contains all the data can be generated. Furthermore, a tar ﬁle is generally
large-sized, so it is better to compress it. The command for compression is gzip cbz.tar, and by doing this we can generate a cbz.tar.gz
ﬁle. The above procedure can be done at once by a command tar cvfz cbz.tgz CBZ, and its opposite is tar xvfz cbz.tgz.
4 Nuclear data deﬁnes probability of interactions between neutron and nucleus. Detail will be described later.
5 Please type cd CBZ/CBG/src.

2

3

Calculations of Takeda benchmark

Here let us calculate a neutron transport problem, Takeda benchmark, with CBZ. 6 The Takeda benchmark includes
four (sub-)problems (or reactor cores), and here we use the simplest one, problem 1. Geometrical speciﬁcation of theproblem 1 is shown in Fig. 2.

Fig. 2: Problem 1 of Takeda benchmark (left: XY-plane, right:XZ(YZ)-plane)
The left part of Fig. 2 is a top view. In this ﬁgure, 1/4 of the whole system is presented. We can ﬁnd that a coreregion is located around the center and a reﬂector region is around the peripheral. Black region in the reﬂector is aposition in which a control rod is inserted. When a control rod is not inserted, this region should be treated as void.
The right part of this ﬁgure is a side view, and you can ﬁnd that a control rod (or void) region has a length along thez-axis, and that a core region is surrounded by a reﬂector region in axial direction also.
Please open a ﬁle main.takeda1.dif 1.cxx in the directory CBGCAL/Takeda.tutorial by an editor.
First you can ﬁnd a comment line with Geometry data around line 20. Below this line, geometry data of problem
is deﬁned. In this problem, Cartesian coordinate is used, and each direction (X, Y and Z) is divided into ﬁve blocks.
As deﬁned by arrays xl, yl and zl, length of each block is set 5 cm, and in calculations, each block is further divided
into ﬁve as deﬁned by arrays xm, ym and zm. Each axis is divided to 25, so the total number of spatial meshes becomes
25 × 25 × 25 = 15, 625 in this example. Material assignment to 5 × 5 × 5 = 125 blocks is deﬁned by an array mat. In
this example, “0” is fuel, “1” is reﬂector, “2” is control rod and “3” is void. In the problem 1 of the Takeda benchmark,two calculations are expected: a control rod insertion case and a control rod withdrawn case. In this example, material“3” is assigned to control rod position, so this example corresponds to the control rod withdrawn case.
Next you can ﬁnd other comment line with the description “Material data” around line 70. Below this line, materialdata, that are neutron-nucleus reaction cross sections, are deﬁned. In this problem, neutron energy is discretized totwo. 7
Finally, using the geometry data and material data deﬁned above, the neutron diﬀusion equation is numericallysolved. As mentioned above, neutron behavior is described by the neutron transport equation. The diﬀusion equation
can be derived by introducing several assumptions to the transport equation. In the neutron transport equation,position, energy and direction dependence of neutrons should be considered, but angular dependence can be neglectedin the diﬀusion equation. This makes diﬀusion equation easier to be solved. Here we use a solver PLOS which cannumerically solve the diﬀusion equation.
First, please type cp main.takeda1.dif 1.cxx main.cxx to copy a ﬁle main.takeda1.dif 1.cxx to main.cxx.
Then type make to compile it. 8
As a result, you can generate a load module ﬁle named a.out. Please conﬁrm this. After the a.out ﬁle is generated,
please type ./a.out and perform calculations. A result similar to the following would be printed on your screen.
6 The Takeda benchmark is a benchmark problem prepared by Prof. Toshikazu Takeda of Osaka University to verify computer programsfor neutron transport. In 1991, experts on neutron transport calculations were invited by OCED/NEA and their own computer codes werecompared with each other using this benchmark problem.
7 Divided energy region is called group, and we can say that the number of energy groups is two in this problem.
8 make compiles a ﬁle main.cxx located in the same directory. Details of procedures are deﬁned in a ﬁle makefile. So when you have aﬁle to be calculated, you have to copy (or rename) it to main.cxx in CBZ.

3

Listing 1: Example of Takeda benchmark calculation
1
2
345678910
11121314151617181920
2122232425262728
2930
31#*** T o t a l Mesh : 15625#*************************************************#* System CBG#*
Power i t e r a t i o n#*
S o l v e r : PLOS#*
A c c e l e r a t i o n : none#*
Forward c a l c u l a t i o n#*
C o n v e r g e n c e c o n d i t i o n ( k e f f ) : 1 e =05#*( flux ): 0.0001#*( source ) : 0.0001#*************************************************#* I t :K eff: Err i n k : Err i n: Err i n#* e r ::: flux: source#*************************************************#0:0.906483
9 . 3 5 2 e =02
2 . 2 3 3 e+00
3 . 7 0 0 e =01#
1:0.917541
1 . 2 2 0 e =02
1 . 0 1 8 e =01
7 . 7 2 3 e =02#2:0.925198
8 . 3 4 4 e =03
6 . 6 3 0 e =02
5 . 3 4 6 e =02#3:0.927640
2 . 6 4 0 e =03
1 . 5 5 2 e =02
1 . 3 7 6 e =02#4:0.927664
2 . 5 6 5 e =05
4 . 6 9 9 e =02
2 . 7 8 7 e =02#5:0.928331
7 . 1 9 8 e =04
4 . 4 2 8 e =03
2 . 9 5 2 e =03#
6:0.927417
9 . 8 4 9 e =04
1 . 4 2 4 e =03
2 . 5 6 3 e =02#7:0.927214
2 . 1 9 2 e =04
4 . 2 0 4 e =03
1 . 6 0 6 e =02#8:0.927261
5 . 1 4 6 e =05
2 . 9 9 1 e =03
7 . 0 6 6 e =03#9:0.927328
7 . 1 2 8 e =05
3 . 1 5 5 e =03
2 . 6 9 0 e =03# 10:0.9273906 . 7 4 4 e =05
6 . 7 1 0 e =04
4 . 4 0 9 e =04# 11:0.927305
9 . 2 2 1 e =05
2 . 3 3 9 e =04
2 . 5 8 0 e =03# 12:0.927291
1 . 4 9 2 e =05
4 . 3 2 4 e =04
1 . 4 4 0 e =03# 13:0.927291
6 . 5 2 0 e =07
8 . 5 5 8 e =04
7 . 0 4 4 e =04# 14:0.927304
1 . 4 1 4 e =05
1 . 4 6 0 e =04
1 . 0 2 2 e =04# 15:
0.927295
9 . 9 2 4 e =06
4 . 2 4 9 e =04
3 . 3 0 9 e =04# 16:
0.927300
5 . 6 8 9 e =06
5 . 6 4 6 e =05
3 . 0 0 9 e =05

A solver PLOS solves the neutron diﬀusion equation with an iterative procedure, so neutron eﬀective multiplication
factor keﬀ and (spatial and energy distributions of) neutron ﬂux are revised at each iteration. By conﬁrming thosequantities are converged, these can be ﬁnal results. More detail, diﬀerence between a result of the present iterationand that of the preceding iteration is compared, and convergence is conﬁrmed if this diﬀerence is smaller than a predetermined criteria. The above example is a process during this iterative procedure, and relative diﬀerences of keﬀ andneutron ﬂux and ﬁssion source are shown. Also you can get the ﬁnal result of keﬀ from this.
Numerical calculations should include various errors and one of the most important ones is discretization error on
space. In the present example, one block with 5 cm length is divided into 5 ×5×5 meshes. We should know that a result
obtained with this model is just an approximated solution, and rigorous solution which is free from the discretization
error can be obtained when we use inﬁnite number of spatial meshes.
Problem 3.1: Perform calculations with ﬁne mesh division: the number of meshes per each plane should be increasedby two or three times. Compare result (keﬀ ) and computation time. You can measure computation time by typing
time ./a.out. Change in the number of spatial meshes can be done by changing values of xm, ym and zm.
This benchmark problem has rigorous numerical solution, which was obtained by solving the neutron transport
equation with the Monte Carlo method, and the rigorous solution is 0.9780 in the control rod withdrawn case. Whatabout your result with PLOS? There might be large diﬀerence between your result and the rigorous one. As alreadymentioned, PLOS yields an approximated solution to the diﬀusion equation, which is approximated form of the transport
equation, so the PLOS solution is just an approximation.
Next let’s use another program to solve the neutron transport equation. Please run main.takeda1.sn.cxx. Notethat this calculations uses 10×10×10 meshes for each block with 5 cm plane. What about your result?
Here we use a solver SNT, which solves the neutron transport equation numerically. By using SNT, you might
get solution close to the rigorous one (0.9780), but you require longer computation time than PLOS. In the ﬁeld ofnumerical calculations, method development to reduce computation time to attain convergence is one of the mostimportant subjects. The same thing can be said in the reactor core analysis ﬁeld, and various methods have beenproposed. Let us use here the diﬀusion synthetic acceleration (DSA) method. This method uses solution of diﬀusion
equation to guess converged solution to transport equation. You can see a command sys.NoDSAAcceleration() aroundline 5 from the end in main.cxx, so please remove this line or comment it out by adding //. By doing this, you canget converged solution with short computation time.
In addition to DSA, we can further accelerate calculation by using the coarse-mesh ﬁnite diﬀerence (CMFD) acceleration. You can use this by changing CalIgen(‘‘none’’); around the end of the program to CalIgen(‘‘cmfd’’);.
Please conﬁrm that computation time is further reduced.

4Problem 3.2: Next let us calculate the control rod insertion case. It can be easily done by changing material
assignment. Control rod reactivity worth can be calculated from two keﬀ s of the control rod insertion and withdrawn
cases, so obtain these reactivity worth by both the diﬀusion and transport calculations. Control rod reactivity
worth can be deﬁned as (1/k1 − 1/k2 ) if neutron multiplication factors before and after the control rod insertion
are k1 and k2 . In the reference Monte Carlo calculations, a control rod reactivity worth of -0.0166∆k/kk ′ wasobtained, so please compare your results with this reference.

Next let’s observe spatial distribution of neutron ﬂux.
Here we use a method ShowNeutronFluxAlongXAxis which outputs neutron ﬂux spatial distribution along the X
axis. You can use this by commenting out a corresponding line in the ﬁle main.takeda1.dif 1.cxx. Arguments ofthis method are as follows: The ﬁrst and second ones are positions in Y and Z directions for which neutron ﬂux spatial
distribution is printed out. If both of them are zero, neutron ﬂux spatial distribution along the X axis about the ﬁrst
meshes on the Y and Z directions. The third argument is an energy group of which neutron ﬂux spatial distribution isprinted out. If you choose -1 for this argument, neutron ﬂuxes of all energy groups are printed. Please be careful thatthe ﬁrst position/energy corresponds to 0 in a C++ array.

Problem 3.3: Please plot spatial distribution of neutron ﬂux of groups 1 and 2 along the X-axis by Gnuplot orExcel. Control rod is inserted and distribution should be chosen over control rod position. In addition, the sameplot should be prepared for the control rod withdrawn case, and compare the distribution before and after thecontrol rod insertion.

Calculation results are printed out on your screen, so you have to copy and paste it. This should be cumbersome, so
you can do calculations by typing ./a.out > output. If you do so, all these results are written in a ﬁle output, andyou can edit this ﬁle to manipulate your results.

5