Resume of PBMR Burnup Analysis
(to be completed)
1. Study by Julian Lebenhaft at MIT (1999)１
In his dissertation at MIT １, Julian Lebenhaft working for the pebble bed
modelling using MCNP4B. One of the analysis is the burnup calculation of the
pebble-bed reactor, and the procedure is applied in PBMR reactor.
The core and fuel element specification of the 268 MWt PBMR used by
Lebenhaft are as follow
Table 1. PBMR Reference Core Specification
Table 2. Reference Specifications for PBMR Fuel Elements
Obara Lab. / Tokyo Institute of Technology | Topan Setiadipura’s note 1
Figure 1. A Vertical View of the PBMR core
To calculate the burnup of the PBMR Lebenhaft utilize a link between
MCNP4B and VSOP94. The fuel management routine in VSOP94 is used to
establish the equilibrium core then nuclide atom densities of the depleted fuel
resulted from this calculation is transfer to MCNP4B as an input. To handle the
data communication between VSOP94 and MCNP4B, Lebenhaft built a program
VSOP94 Model and Analysis
The two-dimensional (R-Z geometry) VSOP94 model of the reference PBMR core is
shown schematically in Figure 2.
Obara Lab. / Tokyo Institute of Technology | Topan Setiadipura’s note 2
Figure 2. VSOP94 Model of the PBR
The model include following detail
- The reflectors and the borated carbon thermal shields
- The reactor pressure vessel, the core barel, and the associated helum gaps are
lumped together in a single annular region.
- The control road are modelled by homogenizing the boron carbide abdorber into
the surrounding graphite material.
- The various channels for the shutdown system and helium coolant are modelled
by reducing the density of the appropriate graphite regions.
- The model of the core comprises of 57 equally sized regions called layers, each of
which contains a mixture of 10 different fuel batches. These batches represent
the number of times that the fuel spheres have been through the core, from the
first pass (batch =1) to last pass (batch = 10.
- The burnup of each batch is calculated individually using a fission-product chain
of 44 isotopes.
- The second member of the chain is a non-saturating fission product chosen by
comparisonwith ORIGEN-JUL-II code predictions, which corresponds to low
absorbing fission products not modelled explicitly.
- The 43 explicit fission products cover 98.02% of the total-fission product
absorptions, while the non-saturating fission product accounts for the remaining
- Fuel depletion is modelled by means of burnup steps.Approximately every 9
Obara Lab. / Tokyo Institute of Technology | Topan Setiadipura’s note 3
days, each layer and its associated batches are moved down to the position of the
next layer. The movement occurs in channels to reproduce the characteristic
flow of pebbles within the core vessel.
- The layers at the bottom of the core discharged, and the corresponding batches
(with the same number of passes through the core) from different channels are
grouped and mixed together. The batches number are the incremented (to
initiate additional pass through the core) and then reinserted into the various
channels. The oldest batches are discharged from the core, while a tenth of the
volume of each top layer is filled with fresh fuel.
- The neutronic calculation is performed using 29 spectral zones, which consist of
user-specified grouping of layers. A similar procedure is used for the reflector
regions where 25 spectrum zones are defined.
- The temperatures of the specified zones are establish using the
thermal-hydraulics module in VSOP94.
MCNP4B Model of the PBMR
The MCNP4B model was based on the 2-D VSOP94 model of the reactor. Key aspect of
the VSOP94 model of the PBMR core and structural components were duplicated,
including the pebble flow channels and layers (Figure 2.).
Figure 2. MCNP4B model of the PBMR (vertical view)
Obara Lab. / Tokyo Institute of Technology | Topan Setiadipura’s note 4
The control rods, and the shutdown and helium coolant in the side reflector
were modelled explicitly with some exeption that the control rods model using the
HTR-10 model (instead of proprietary PBMR design), also the oval shutdown channels
which taken from Chinese reactor.
A simple BCC lattice without exclusion zones was used to represent the
randomly packed pebble bed. This lattice was subdivided into layers exactly as in the
VSOP94 model of the core. The spheres in each layer consisted of the homogenized mix
of carbon, silicon, heavy metal dan fission-product nuclides prepared with MCARDS
from the corresponding VSOP94 layer-averaged atom densities.
Link Between VSOP94 and MCNP4B
An approximate but direct approach was used by Lebenhaft where the nuclide densities
obtained from a whole-core VSOP94 calculation are used directly in an MCNP4B model
of the reactor.
As explained before, an MCNP4B model of the PBMR equilibrium core was
developed, which uses a regular lattice of spheres with homogenized fuel and graphite
interiors. This lattice is subdivided into annular regions that reproduced exactly the
core regions of the VSOP94 model. The fuel composition used in these regions are
batch-aveeraged VSOP94 muclide densities, which are adjusted to allow for the effect of
fuel homogenization on resonance absorption in 238U and 240Pu.
In the VSOP94 code,subroutine VORSHU is called before fuel shuffling takes
place and work the following jobs
- Reads instructions for the next burnup
- Calculates average batch properties
- Account for isotope decay during the reshuffling.
this VORSHU subroutine was midified to accommodate new jobs as follow
- Calculate the total atom numbers and write the result to FORT.98
- Calculate the reaction rates for each isotope per layer and write the result to
- Additional information written to FORT.98 file are
Total number of all nuclides in the given layer
The volume of the layer
The packing factor.
A FORTRAN program MCARDS reads these data and prepares an MCNP4B
Obara Lab. / Tokyo Institute of Technology | Topan Setiadipura’s note 5
material card for each layer. Because a one-to-one mapping does not exist between the
nuclides available in the VSOP94 and MCNP4B libraries, is 10B added such that the
total absorption reaction rate in each layer is preserved.
With the above hybrid method VSOP94 performs the diffusion calculations on
homogenized regions of the core, with self-shielding factors used to correct for double
heterogeneity and a diffusion correction to allow for neutron streaming in core voids.
Neutron streaming is handled in MCNP4B by modelling explicitly the lattice of spheres.
Preliminary analysis of the PBMR using the above MCNP4B/VSOP94
calculation method yields an effective multiplication constant at 1.00873 ± 0.00096. The
nuclide densities calculated by VSOP94 for each layer is not adjusted, this adjustment
is needed to compensate for the approximate MCNP4B representation of fuel-sphere
interiors. The effect of this homogenization was shown to decrease the reactivity of the
core ~4%dk/k by using an identical MCNP4B model in which the core layers were fully
homogenized. Two-dimensional VSOP analysis of the PBMR core has established a total
temperature reactivity shown in the table 3. Which would reduce the reactivity of the
core by ~4%dk/k assuming an average fuel temperature of 950oC. The MCNP4B
calculations reported were performed at room temperature, because cross-section
libraries at typical HTGR operating temperatures were not available.
The method described in this chapter for obtaining the MCNP input for the
composition of fuel spheres in a pebble-bed core with burnup islimited in its accuracy,
because of the need to account for the effects of homogenization on resonance absorption.
To improve the calculation method, Lebenhaft proposed the following procedure for
MCNP-based code system.
Obara Lab. / Tokyo Institute of Technology | Topan Setiadipura’s note 6
Figure 3. Proposed MCNP-bsed Code System for Pebble-Bed Reactors with
2. Study by Brian Boer et.al at TU Delft (2009)2,3
In the paper published 20092, Boer et.at evaluating several pebble loading
schemes.The evaluation is based on inner and outer graphite reflector which
included in some of current pebble-bed reactor designs with increased thermal
power (>300 MW) such as the designs considered in the HTR-PM (China) and 400
MW PBMR (South Africa).
The above annular pebble-bed design can exhibit peaks near this reflectors
caused by the local abundance of thermal neutrons. Furthermore, the tall core
geometry, adopted for thermal-hydraulic reasons, causes a large difference in the
burnup level between top and bottom of the core, resulting in an axial power peak at
Obara Lab. / Tokyo Institute of Technology | Topan Setiadipura’s note 7
the top. At the end, this power peaks can results in hig fuel temperatures during
normal operating and accident conditions.Thus, optimized loading pattern is needed
to improve the power profile which can lead to a reduction of the fuel temperature
for a fixed helium outlet temperature.
To evaluate several loading pattern effect on the reactor parameters, Boer et.al
developed a new calculation system to be expalained here and applied the procedure
to 400MWt PBMR.
Reference Design of 400 MWt PBMR
Reference design of the 400 MWt PBMR used for the evaluation is as follow
Figure 2.1 General design and operating characteristic of the PBMR 400 MWt.
Koster4 give a schematic picture of the core design as follow
Figure 2.2 Core configuration
In the calculation it is assumed that the equilibrium core is already achieved so that the
Obara Lab. / Tokyo Institute of Technology | Topan Setiadipura’s note 8
neutron flux profile does not change in time. Also assuming that the pebble is irradiated
wih a fixed flux level during a certain time interval, when it moves from one point in the
core to the next.
The general equation for fuel depletion, including fuel movement is as follows
ji aj ji j N j ai j N i
dt j 1
= Material derivative of nuclide i
Ni = Atomic concentration of atom i
= Neutron flux
ji = Probability that neutron interaction with nuclide j will yield nuclide i
ji = Probability that the decay of nuclide j will yield nuclide i
j = Decay constant for nuclide j
aj = Absorption cross section of nuclide j
The outline of the calculation method id depicted in figure 2.2 as follow
Obara Lab. / Tokyo Institute of Technology | Topan Setiadipura’s note 9
Figure 2.2 Flow scheme of calculation procedure for the calculation of the power profile
in the core.
The above procedure involving several codes. The codes are used iteratively unti
convergence is reached on both the inner and the outer iteration. The criterion for the
inner loop is the change in the neutron flux, while convergence in outer loop is
determined by the deviation from a preset kdes value. The final burnup level of the
pebbles, Bufin, id modified to meet this criterion.The total residence time T that the
pebbles reamin in the core is modified in turn to reach this burnup level.
The inner loop consist of the following step :
1. Depletion calculation.
- Provided the flux profile in the core is known the fuel in the pebble during its
lifetime can be calculated.
- The core is divided into several discrete radial and axial zones
- Depletion calculation is performed using the ORIGEN module from SCALE-5.
- Between burnup steps the cross sections are updated using successive 1D
Obara Lab. / Tokyo Institute of Technology | Topan Setiadipura’s note 10
transport calculations for the TRISO coated particle and pebble geometry.
- 172 energy group XMAS library is used based on the JEFF3.1 libarary.
- Results of this step is the nuclide distribution over the core of each burnup
2. Neutron cross section calculation
- Zone averaged nuclide concentrations N(r,z) are used for generating neutron
cross section for the entire reactor.
- Double heterogeneity of the fuel is taken into account by performing a cell
weighting of a fuel particle in moderator material using a Dancoff factor
calculated according to the method developed by Bende5.
- Resonance treatmen of the unresolved resonance is performed by Bondanrenko
method, While the resolved method using the Nordheim Integral Method.
- Then followed by a 1-D transport calculation of a pebble geometry including the
surrounding helium. Radial and axial (1-D) transport calculations are used to
calculate zone weighted cross section for the entire reactor.
- Result of this step is a two dimensional cross section library (Σ(r,z)).
3. Neutron flux calculation
- The above neutron cross sections are used to calculate the 2-D multigroup flux
profile (r,z) with the neutron diffusion code DALTON.
- Average fluxes are generated for each depletion zone and scaled to the desired
- K-eff of the reactor is calculated
- Using the nuclide distribution over the core and the cross sections are generated
for several temperatures of the fuel and moderator, the temperature dependent
cross section library is generated.
- The library can be used in coupled neutronic and thermal hydraulic calculations
using the DALTON-THERMIX code system. The results is a temperature
corrected power profile.
1. Lebenhaft JR, 1999, “MCNP4B Modelling of Pebble Bed Reactor”, MIT Dissertation.
2. Boer B, Kloosterman J.L, lathouwers D, van Der Hagen T.H.J.J, 2009, “In-core fuel
management optimization of pebble-bed reactors”, Annals of Nuclear Energy 36
(2009) 1049 – 1058
Obara Lab. / Tokyo Institute of Technology | Topan Setiadipura’s note 11
3. Boer B, 2009, “Optimized Core Design and Fuel Management of a Pebble-Bed Type
Nuclear Reactor”, TU Delft Dissertation.
4. Koster A, Matzner H.D,NIcholsi D.R, 2003, “PBMR Design for the Future”, Nuclear
Energy and Design 222 (2003) 231 – 245
5. Bende E.E,1999,”Pebble Bed Burning in a Pebble Bed Type High Temperature
Nuclear Reactor”, TU Delft Dissertation.
Obara Lab. / Tokyo Institute of Technology | Topan Setiadipura’s note 12