Modeling of realistic pebble bed reactor geometries using the Serpent Monte Carlo code

Modeling of realistic pebble bed reactor geometries using the Serpent Monte Carlo code

Annals of Nuclear Energy 77 (2015) 223–230 Contents lists available at ScienceDirect Annals of Nuclear Energy journal homepage:

2MB Sizes 0 Downloads 63 Views

Annals of Nuclear Energy 77 (2015) 223–230

Contents lists available at ScienceDirect

Annals of Nuclear Energy journal homepage:

Modeling of realistic pebble bed reactor geometries using the Serpent Monte Carlo code Ville Rintala a,⇑, Heikki Suikkanen a, Jaakko Leppänen b, Riitta Kyrki-Rajamäki a a b

LUT Energy, Lappeenranta University of Technology, P.O. Box 20, FI-53851 Lappeenranta, Finland VTT Technical Research Centre of Finland, P.O. Box 1000, FI-02044 VTT, Finland

a r t i c l e

i n f o

Article history: Received 18 July 2014 Received in revised form 29 October 2014 Accepted 17 November 2014

Keywords: High temperature reactor Pebble bed reactor Reactor physics Monte Carlo Discrete element method

a b s t r a c t This paper documents the models available in Serpent for high temperature reactor (HTR) calculations. It is supplemented by a calculation example of ASTRA critical pebble bed experiments. In the pebble bed reactor modeling, different methods have been used to model the double heterogeneity problem occurring in pebble bed reactor calculations. A solution was sought to avoid unphysical simplifications in the pebble bed modeling and the stochastic geometry modeling features available in the Monte Carlo code Serpent were applied for exact placement of pebbles and fuel particles. Randomly packed pebble beds were produced in discrete element method (DEM) simulations and fuel particles were positioned randomly inside the pebbles. Pebbles and particles are located using a Cartesian search mesh, which provides necessary computational efficiency. Serpent uses Woodcock delta-tracking which provides efficient neutron tracking in the complicated geometries. This detailed pebble bed modeling approach was tested by calculating the ASTRA criticality benchmark experiment done at the Kurchatov Institute in 2004. The calculation results are in line with the sample calculations provided with the benchmark documentation. The material library selected for the calculations has a major effect on the results. The difference in graphite absorption cross section is considered the cause of this result. The model added in Serpent is very efficient with a calculation time slightly higher than with a regular lattice approximation. It is demonstrated that Serpent can be used for pebble bed reactor calculations with minimal geometric approximations as it allows exact pebble bed modeling with randomly positioned fuel particles and locations of pebbles produced by DEM. Due to the added stochastic geometry features and efficient neutron tracking, it is concluded that Serpent is well suited for the calculations of HTRs. Ó 2014 Elsevier Ltd. All rights reserved.

1. Introduction The ability to define reactor geometries with a high level of detail make Monte Carlo method based reactor physics especially attractive for modeling pebble bed reactors (PBR), where the fuel is enclosed inside millimetre scale coated particles and dispersed inside graphite spheres. The pebbles themselves are randomly packed inside the reactor core. Traditionally Monte Carlo models of pebble beds have been simplified for computational efficiency by having the pebbles as well as the fuel particles in regular lattices and sometimes experimental configurations have been performed with regularly stacked pebbles to allow easier modeling and direct comparison between calculations and experiments (Difilippo, 2003). A real PBR however, does not consist of regularly placed pebbles. In Serpent the preferred approach is to explicitly model the stochastic positions of both the pebbles and the fuel particles inside ⇑ Corresponding author. Tel.: +358 40 152 9437. E-mail address: [email protected]fi (V. Rintala). 0306-4549/Ó 2014 Elsevier Ltd. All rights reserved.

them. A good computational performance is obtained by using the Woodcock delta tracking method in conjunction with the explicit stochastic geometry models implemented in Serpent. For the generation of realistic dense pebble beds, the discrete element method (DEM) is suggested. DEM is a well established computational method for modeling particulate systems that has been used for pebble bed reactor related analyses previously, e.g., by Suikkanen et al. (2014) in pebble packing studies. In this paper, a methodology for modeling pebble bed reactors with Serpent is explained together with the use of DEM for constructing realistic pebble beds. To demonstrate the methodology, pebble bed criticality experiments performed in Russia between 2003 and 2004 with the ASTRA facility are calculated. The ASTRA benchmark provides a challenging modeling task with its geometrical peculiarities and stochastic pebble cores. The ASTRA benchmark configurations have previously been calculated with Serpent and documented in Suikkanen et al. (2010). In the work presented in this paper, however, major improvements have been made to the calculation model. These include inclusion


V. Rintala et al. / Annals of Nuclear Energy 77 (2015) 223–230

of the geometric details of the bottom reflector and completely new DEM simulations to obtain stochastic pebble bed configurations that have the core heights within the measurement tolerances reported for the experimental configurations. Also a more complete description of the methodology for modeling the double heterogeneous stochastic pebble bed geometries in Serpent is given. 2. Methodology in the Serpent code One of the advantages of using the continuous-energy Monte Carlo method for particle transport problems is the ability to model complicated geometries without major approximations. Complex structures in reactor physics applications are typically related to fuel assembly geometries in which the fuel pins are arranged in a regular lattice structure. The use of regular lattices not only simplifies the model construction for the user, but also considerably improves the efficiency of the tracking routine. When a neutron path is started from the coolant region between fuel pins, the tracking routine needs to determine the next boundary crossing within the line-of-sight, which is found by calculating the distances to the pin surface and the unit cell boundaries. If the same geometry structure is modeled explicitly without the regular lattice, the tracking routine would have to loop over all the fuel pins within the assembly or reactor core in order to find the next boundary crossing. Since the number of fuel pins in full-core geometries increases to several thousand, the tracking routine would become a major bottleneck for the entire Monte Carlo simulation. Even though regular lattice structures work extremely well for most fuel and reactor types, they are unable to capture the randomness effects of single- and double-heterogeneous stochastic geometries. There exist methods for introducing randomness in regular lattices (Brown and Martin, 2004), but the best way to model full-scale pebble-bed type high temperature reactors (HTR) is to put each fuel particle and pebble explicitly in its place. The explicit description makes it possible to not only capture the randomness effects, but to also avoid particle clipping at the outer boundary – the two factors that must be taken into account in order to avoid systematic errors in HTR modeling. The problem with the explicit approach is that the efficiency of the Monte Carlo tracking routine is completely lost. The number of fuel particles inside a typical pebble is over 10,000, and the number of pebbles in the core ranges from 27,000 in a test reactor (HTR-10) (Wu et al., 2002) to 675,000 in a full-sized reactor (THTR) (Melese and Katz, 1984). If the tracking routine needs to calculate the distance to every surface each time a path length is sampled, the cost in CPU time becomes prohibitively high and the method too impractical to use. 2.1. Explicit stochastic geometry model Division of complicated geometry structures into smaller units is a well-known method for improving the efficiency of the Monte Carlo tracking routine, and this is exactly what is done automatically when a conventional reactor geometry is modeled as a regular lattice of assemblies and pins. The same approach was taken as the basis for the explicit stochastic geometry model developed in the Serpent code. The method has been available in Serpent 1 since the public release of version 1.1.0 in March 2009, but its use has not been particularly well documented.1 The same method is also implemented in the development version of the code, Serpent 2. 1 It should be noted that Serpent 1 also has a second, implicit particle fuel model (Leppänen, 2007), which is based on the random sampling of fuel particles along the neutron path. The implicit model, however, fails to account for the finite size of the fuel particles, which results in significant errors when the packing fraction is high. The method is therefore not recommended for use when modeling HTR geometries.

Fig. 1. Section of modeled pebble bed with both pebble and particle level search meshes of Serpent visible.

The explicit model reads the coordinates of spherical objects, such as fuel particles or pebbles, from a separate input file and creates a three-dimensional lattice-like structure that can be used as any other universe in the Serpent geometry. Instead of fitting the objects in a regular pattern, the code creates a Cartesian search mesh over the distribution. Every mesh cell has a list of pointers to spheres that are located inside the cell or intersect its boundaries. Tracking is carried out from one mesh cell to the next, similar to a conventional pin-cell lattice, and the geometry routine needs to check boundary crossings only with the spheres that have pointers assigned to that cell. The number of boundary checks can be reduced by several orders of magnitude, which results in a considerable increase in calculation speed. Example of this Cartesian search mesh can be seen in Fig. 1. The model is not based on any approximations, and it can be used for modeling both particle distributions inside the fuel pebbles and pebble distribution in the core.2 Generation of random particle distributions for the input files can be carried out using a separate disperser routine in the Serpent code. Methods for generating realistic pebble distributions are discussed in Section 3. In addition to input, the same structures can be used for reaction rate tallies. If requested by the user, pebble-wise power distribution is calculated automatically and printed in a separate output file together with the explicit coordinates of each pebble. 2.2. Woodcock delta-tracking method Another problem related to HTR modeling with continuousenergy Monte Carlo codes is the fact that the neutron mean-freepath is long compared to the spatial dimensions of the system. In a regular lattice model the neutron may cross several cell boundaries before reaching the collision site, and the tracking routine needs to stop the path at each boundary crossing. The same problem occurs with the search mesh of the explicit stochastic 2 Even though the model is not based on any approximations, the scale of most HTR geometries necessitates the use of the same random particle distribution for multiple spheres or compacts. This is not considered a major problem, due to the large number of particles, and no biases have so far been observed in any test calculations.

V. Rintala et al. / Annals of Nuclear Energy 77 (2015) 223–230

model. The overall result is that a significant fraction of CPU time is spent merely calculating the distances to the next cell boundary. To avoid this loss of efficiency, the tracking routine in Serpent is based on a combination of conventional surface-tracking and the Woodcock delta-tracking method (Leppänen, 2010). Delta-tracking is essentially a rejection sampling technique that is used to transport the neutron from one collision site to the next without stopping the track at surface crossings. In LWR calculations this technique can speed up the tracking routine by a factor of two, but in HTR particle fuel geometries the factor can be assumed to be substantially higher. For the explicit stochastic geometry model the use of deltatracking means that the search mesh can be refined indefinitely without loss of efficiency due to shortened track lengths. Experience has shown that the best result is obtained when the cell size is set comparable to particle diameter, which means that the number of particles per cell is minimized. The running time then becomes comparable to a regular-lattice model with the same packing fraction, and the randomness of the distribution can be accounted for without penalty in efficiency. 3. Construction of realistic pebble beds The volume fraction of a mechanically stable random packing of mono-sized spheres varies from 0.55 to 0.64 (see, e.g. Onoda and Liniger (1990) and Song et al. (2008)). The average volume fraction in the bulk region (away from walls) of a pebble bed core should lie between these values but in reality it is more likely between 0.609 and 0.625 when the spheres are poured into a container (Dullien, 1979). Higher values can be obtained when the packing is vibrated. Because of the high volume fraction, it is impossible to generate random pebble bed configurations with a simple random number generation algorithm in the same way as particle distributions inside the pebbles as described in Section 2.1. Either a more sophisticated packing algorithm or a method that models the physical core loading process itself is needed. Several sphere packing algorithms exist; however they are mostly based only on geometric parameters and thus do not resolve the actual contact mechanics between the spheres. See, e.g., work by Visscher and Bolsterli (1972) and Jodrey and Tory (1985). Although such algorithms can be used for generating random pebble distributions with a minimal computational effort, great care must be taken to verify that the pebble beds actually have physically realistic packing structures. These algorithms have a limited applicability and cannot be used for further dynamic analyses, e.g. simulations of pebble flow or behavior of the pebble bed during an earthquake. 3.1. Discrete element method A physically realistic modeling approach using discrete element method (DEM) introduced by Cundall and Strack (1979) is suggested for the construction of pebble beds to be used in Serpent calculations. DEM is closely related to molecular dynamics methods as it uses some similar physical models and algorithms. In DEM, the trajectories of discrete particles, typically spheres, are calculated as they interact with each other and domain boundaries and are affected by external forces such as gravity. Particles are regarded as soft spheres that interact when they overlap. Normal and tangential contact forces between a particle pair are calculated with a suitable force model that is formulated as a function of their overlap distance that represents the deformation between the particles. Several contact force models exist, see, e.g., Di Renzo and Di Maio (2004), and they usually consist of spring and dashpot mechanical components. Both translational and rotational


movement of the particles can be considered and the force model coefficients can be formulated based on the material parameters of the particles. At each DEM simulation time-step, the net force on each particle is calculated and the Newton’s equations of motion are solved. Due to usually quite high stiffness of the particles, the repulsive force in a contact grows steeply. This requires a rather small time-step size. A DEM code also needs to update a list of the neighbors of each particle. These features make DEM a computationally expensive method. Nevertheless, the popularity of DEM has been constantly increasing due to the increase in computer power and it is used in modeling of discontinuous materials in many fields. In nuclear engineering, DEM can be used to construct realistic pebble beds inside any core geometry. It is also capable of modeling dynamic behavior such as pebble flow as demonstrated, e.g., by Rycroft et al. (2006). Lots of useful data can be extracted from a DEM simulation including, e.g., the forces exerted on reflectors by the pebbles. Although computationally demanding, the simulation of pebble beds with DEM is feasible as the typical number of pebbles is relatively small compared to more conventional application areas of DEM where typically systems consisting of millions of particles are simulated. Several commercial and free DEM codes are available with varying number of implemented force models and other features, such as capability to read in CAD geometries, interfaces and coupling with fluid flow solvers, and parallel processing capabilities. 4. Calculation model of ASTRA critical experiments The ability of Serpent to model PBRs is demonstrated by calculating the five critical pebble bed configurations that were assembled in experiments with the ASTRA facility in the Russian Research Center (RRC) Kurchatov Institute between 2003 and 2004 and documented by Ponomarev-Stepnoi et al. (2008). With stochastic distribution of pebbles, the ASTRA experiments provide an excellent benchmark to demonstrate the pebble bed models available in Serpent. A brief description of the facility and the experiments is given below. A more detailed description can be found in the original benchmark documentation (PonomarevStepnoi et al., 2008) or in Garin et al. (2010). The short introduction to the experiments is followed by a description of the detailed calculation model built for Serpent including the construction of the pebble beds using DEM. 4.1. Description of the facility and the experiments The ASTRA critical configurations were assembled inside a cylindrical steel vessel with an outer diameter of 3.82 m and a height of 4.63 m. Side, bottom and center reflectors were assembled inside the vessel from graphite blocks with a square cross section. Each block had a small recess at the bottom and a corresponding projection on the top to help alignment when piled on top of each other. The blocks also had a cylindrical central hole, which was filled with a graphite plug in positions where they were not used as channels for control rods (CR), detectors or other instrumentation. The side reflector was shaped in the form of an octagon (side lengths of 175.56 cm in horizontal and vertical directions and 178.19 cm in diagonal directions). This was done by cutting some of the graphite blocks diagonally along the vertical axis. The holes in these blocks were not plugged, which resulted in small recesses in the core cavity. The center reflector was shaped as an octahedral prism (side lengths of 99.54 cm in horizontal and vertical directions and 105.73 cm in diagonal directions) by cutting some of the blocks in the same way as in the side reflector. Unlike in the side reflector, these blocks were plugged. Other center reflector blocks did not


V. Rintala et al. / Annals of Nuclear Energy 77 (2015) 223–230

with natural boron carbide. In addition, a single manual control rod (MR) was used that contained only air inside an aluminium tube. The MR had its position in the side reflector. The five experiments had a different number of fuel pebbles, resulting in different core heights. The poison rods were also at different positions. In the last experiment, a top reflector consisting of a stochastic bed of graphite spheres similar to the fuel spheres was placed on top of the core with a 2 mm thick aluminium sheet in between. The main parameters of the experiments are given in Table 1. 4.2. Pebble bed configurations

Fig. 2. Approximate schematic of the ASTRA annular core assembly showing the various instrumentation and reactivity control positions.

Table 1 Main parameters of the ASTRA configurations. Core

Fuel/moderator pebbles

MR position (cm)

CR 5 position (cm)



1 2 3 4 5

16,897 20,287 27,671 30,432 30,432/9512

178.8 160.5 225.1 Out Out

Out Out Out 184.6 93.0

Out In In In In

Out Out In In In

have plugs and some of the channels were used for leave-in-place poison rods (LIPR). An additional square cross section channel was made for CR 7 poison rod and several smaller holes were drilled in some blocks for detectors or profiling absorber elements. The bottom reflector consisted of a single layer of blocks (height of 40.00 cm), which were all plugged. Together the reflectors formed an octagonal annular cavity (height of 420 cm and thickness of 38.01 cm in horizontal and vertical directions and 36.23 cm in diagonal directions) where the pebble bed core could be assembled. A descriptive schematic of the facility showing the locations of the various control rods and instrumentation is given in Fig. 2. The cavity was filled with a variable number of spherical fuel elements depending on the experiment. The fuel elements were of typical design with an approximate outer diameter of 60 mm. The coated fuel particles were randomly dispersed in graphite inside the fuel zone with a diameter of 50 mm. The fuel particles had a uranium dioxide kernel surrounded by a layer of porous pyrolytic carbon, and two layers of dense pyrolytic carbon with a layer of silicon carbide in between. Each fuel element contained 2.44 g of uranium. The enrichment in 235U was about 21%. A channel for CR 6 was placed inside the cavity. Additionally six detector channels were placed vertically along the core mid-plane and eight detector channels were placed near the side reflector boundary inside the half cylinder recesses. The control rods, safety rods and leave-in-place poison rods used in the experiments were of the same design. Each rod was an assembly of steel tubes filled

An open source DEM software LIGGGHTS (LAMMPS Improved for General Granular and Granular Heat Transfer Simulations, (Kloss et al., 2012) is used for the construction of the pebble beds in this work. LIGGGHTS has been built on top of the molecular dynamics simulation software LAMMPS (Large-scale Atomic/Molecular Massively Parallel Simulator, (Plimpton, 1995). LIGGGHTS is actively developed and is suitable for modeling of pebble beds. It has parallel computing capabilities via MPI (Message Passing Interface) and complex wall geometries can be imported in STL (Stereolithography) format. For the ASTRA packing simulations, a detailed representation of the walls surrounding the pebbles was built using external modeling tools, in this case the computational fluid dynamics pre-processing software Gambit by ANSYS Inc. and exported as an STL file. The wall geometry used in the DEM simulation is shown in Fig. 3. In addition to the boundaries defined by the center, side and bottom reflectors, the model includes the in-core control rod channel and all detector channels located inside the core. The channels are defined using the cylinder wall primitive available in LIGGGHTS. Projections on top of the bottom reflector graphite blocks are included in the model for two reasons. First, they have an effect on the pebble packing structure in the bottom region as their height is almost half of a pebble radius. Second, because the core height in the experiments was calculated as the distance from the top of these projections to the top of the pebble bed. No detailed description of the core loading process in the experiments is given by Ponomarev-Stepnoi et al. (2008). It is assumed that the core was filled with pebbles using some kind of feed tubes and that the surface of the pebble bed was flattened to get an even core height throughout the cavity. To mimic a real core loading process, the pebbles in the packing simulation were dropped from four inlets located in the north-east, south-east, south-west and north-west parts of the cavity. Ten pebbles per inlet were created at a time and let to fall due to gravity. In preliminary simulations it was found that this packing method produces slightly too high pebble beds, i.e. too dilute packings. To get pebble beds closer to those reported in the experiments, the domain walls were subjected to a high-frequency oscillating motion after the compaction of the initial packing so that the vibration induces additional densification in the pebble bed. Thus the simulation has four phases. In the first phase, a total number of 36,000 pebbles are packed inside the cavity by dropping them from the inlets. In the second phase, the pebble bed is vibrated for a suitable amount of time. In the third phase, a number of top pebbles is removed corresponding to each experimental case, leading to an even pebble bed surface and the correct amount of pebbles for the particular core configuration. In the fourth phase, the simulation is continued to let the modified packed bed relax to a stable state. In the last ASTRA core configuration, an additional wall is defined on top of the highest fuel pebble of the fourth core configuration and a layer of moderator pebbles is formed on top of the wall in the same way as the fuel pebble cores are formed.

V. Rintala et al. / Annals of Nuclear Energy 77 (2015) 223–230

Fig. 3. ASTRA wall geometry used in the DEM packing simulation. Pebbles are packed inside the annulus between the center and side reflectors. The in-core control rod channel (black circle), detector channels (small black dots) and the 14 mm high ring-shaped projections on top of the bottom reflector blocks (grey circles) are included in the model.

Table 2 Parameters used for the pebbles in the DEM packing simulation. Parameter



Fuel and moderator pebble diameter Fuel pebble density Moderator pebble density Young’s modulus Poisson’s ratio Restitution coefficient Friction coefficient

5:985  0:002 1908 1681 10 0.3 0.9 0.3

cm kg/m3 kg/m3 GPa – – –

Realistic pebble densities and diameters are used given by Ponomarev-Stepnoi et al. (2008). Suggestions given in literature for other mechanical properties, e.g. the friction coefficient by Xiaowei et al. (2010), are used. The values used are given in Table 2. A Gaussian distribution of pebble diameters using the mean value and the standard deviation given by Ponomarev-Stepnoi et al. (2008) is used. The small variation in the pebble diameters is not considered to have any real effect in either the DEM or the Serpent calculations, but is included to emphasize how each pebble can be made unique. The density values are used for calculating pebble masses. In the case of fuel pebbles the value is for the homogenized mixture of fuel particles and graphite. 4.3. Serpent model A detailed full core reactor model was built for Serpent with a stochastic pebble bed and randomly positioned fuel particles inside the pebbles. With an accurate model of the pebble bed itself, no major approximations were needed to model the benchmark cases. A description of the ASTRA facility and the experiments was given above and a more detailed description is available in the benchmark documentation. In the following sections, only


the simplifications to the actual case and differences to the benchmark specifications in our model are given. Pebble beds were modeled with the explicit method described in Section 2.1 and from the coordinates obtained from the LIGGGHTS simulations. A feature of Serpent to allow pebbles with different radiuses to be defined was used to model a small variation in the sizes of the pebbles. This was already taken into account during DEM simulations. This variation of radiuses changed the amount of graphite only in the shell covering the fueled zone of the pebbles. In the DEM simulations the pebbles overlap each other slightly, which leads to a small graphite loss in the reactor physics calculations. The amount of lost graphite due to this modeling artifact was calculated to be insignificant when compared to the total amount of graphite in the pebbles. All details visible in Fig. 3 in the core cavity modeled in LIGGGHTS were also modeled in Serpent. This included the graphite projections on top of the bottom reflector graphite blocks, instrumentation tubes inside the core, and the tube of the control rod CR6. The bottom attachment of the control rod tube CR6 on the bottom of the core cavity was not modeled. Slanted surfaces of the ring shaped projections on top of the graphite blocks were straightened but the amount of graphite was kept unchanged. This applies to all projections which were modeled in a detailed way throughout the geometry. Fuel particles inside the pebbles were modeled with the average dimensions reported in the benchmark documentation, all particles being similar throughout the geometry. The particles were reported to contain a certain mass fraction of impurities without more precise information about their positions in the particles. Therefore the impurity mass fractions were assumed to be constant in the uranium kernel and all the different coating layers. The particle locations were generated with the disperse feature and positioned with the explicit stochastic geometry model in Serpent. Ten different particle distributions were generated and these distributions were randomly chosen for the pebbles. Fig. 1 shows a cross section of several fuel pebbles with randomly distributed particles and the search meshes used by the geometry routine of Serpent. In general, all reflector graphite blocks were modeled as solid graphite and air gaps between the blocks were taken into account by modifying the graphite density. All graphite blocks had an axial hole going through them and some of them were filled with graphite plugs. Empty holes were modeled as such, but the holes with filler plugs were treated as solid graphite. The internal reflector was standing on a support structure made from aluminium and was modeled as reported in the documentation. Instrumentation channels and experimental grooves inside the center reflector were included in the model but the top projections of the graphite blocks were leveled. Graphite blocks were also leveled on top of the side reflector. On the bottom reflector the recesses were only modeled by changing the overall graphite density. Experimental grooves in the side reflector were not modeled. Anything above top surfaces of internal and side reflector were not included in the model. This includes the graphite column lids made from aluminium, control rod tubes and other structures on top of these surfaces. Neutron detectors, cabling and other instrumentation related equipment were not modeled. In the model, neutrons leak from the top when the level 500 cm from the base of the bottom reflector is reached. In the top part of the reactor core details were modeled accurately to this level. The base of the bottom reflector is the limit in the opposite direction. Sideways the model was limited to the inner surface of the stainless steel vessel. A horizontal cross section of the Serpent model in benchmark case 5 is presented in Fig. 4 and vertical cross section in Fig. 5. The material libraries used in the calculations were JEFF-3.1.1 at 300 K, JENDL-4 at 300 K and TENDL-2012 at 293 K. Thermal scattering data were used for graphite and the used library was JEFF-3.1 at 294 K. Impurities in graphite were modeled by adding one ppm by


V. Rintala et al. / Annals of Nuclear Energy 77 (2015) 223–230 Table 3 Comparison of the core heights in the experiments and in the calculation model. Core

Core height in experiment (cm)

Core height in simulation (cm)

Difference (%)

1 2 3 4 5

179:36  0:53 214:14  0:53 291:59  0:48 320:05  0:41 421:58  0:19

179.52 214.04 291.57 319.92 421.46

0.09 0.05 0.01 0.04 0.03

Fig. 4. Horizontal cross section of the Serpent model of the ASTRA benchmark case 5.

Fig. 6. Pebble bed corresponding to ASTRA core configuration 1 with 16,897 pebbles constructed using DEM.

5. Results and discussion

Fig. 5. Axial cross section of the Serpent model of the ASTRA benchmark case 4.

weight of natural boron as instructed in the benchmark documentation. Material data for stainless steel and aluminum alloys were taken from Russian GOST standards. No specific information about neutron absorbing impurities were available for these materials.

The stochastic pebble beds obtained from the DEM simulations were evaluated and compared to the experimental configurations. In the experiments, the core height was determined as an average of height measurements at eight positions measured with a specific measurement device (Ponomarev-Stepnoi et al., 2008). A similar method was used to determine the heights of the numerical pebble beds. The heights of the pebble beds produced using DEM and the corresponding experimental core heights are given in Table 3. It can be seen that the core height of each configuration is inside the given measurement uncertainty of the experiments. An example of a pebble bed resulting from a DEM simulation is shown in Fig. 6. The total number of neutron histories simulated in each case was 252 million with 500 actual calculation cycles. The results of the first ten cycles were disregarded and these neutron histories are not counted for total. The experimental results are given in Table 4 together with the sample calculation provided with the benchmark and our results calculated with Serpent. Results were calculated with a computer cluster using Intel Xeon E5-2660 processors and the total CPU time used was 300 h or less depending on the benchmark case. Number of parallel tasks used were 64 or less. The explicit stochastic geometry model in Serpent performed as expected and allowed accurate modeling of the ASTRA criticality facility. This calculation case and other similar cases have shown


V. Rintala et al. / Annals of Nuclear Energy 77 (2015) 223–230 Table 4 Multiplication factors for different ASTRA configurations. Core

keff  1r (experiment)

keff  1r (sample results) MCU-REA1 DLC/MCUDAT-2.2

keff  1r (simulation) Serpent 1.1.19 JEFF-3.1.1

keff  1r (simulation) Serpent 1.1.19 TENDL-2012

keff  1r (simulation) Serpent 1.1.19 JENDL-4

1 2 3 4 5

1:0000  0:0036 1:0000  0:0036 1:0000  0:0036 1:0000  0:0036 1:0000  0:0036

0:9912  0:0005 0:9936  0:0005 0:9977  0:0005 0:9989  0:0005 1:0006  0:0005

1:00805  0:00008 1:01038  0:00008 1:01216  0:00008 1:01305  0:00008 1:01403  0:00008

0:99922  0:00008 1:00256  0:00008 1:00537  0:00008 1:00633  0:00008 1:00772  0:00008

0:98963  0:00008 0:99292  0:00008 0:99520  0:00008 0:99620  0:00008 0:99726  0:00008

that typically calculation times have been 10% or less longer compared to cases where the pebbles and particles have been in a regular lattice. Overall, the penalty caused by the accurate modeling approach seems to be minimal. The number of disregarded cycles in the beginning of each calculation were set to 10 and fission source entropies were checked to ensure the converged fission source. Depending on the amount of modeled detail there is bias between the modeled and measured results. In the benchmark documentation the effect of the different factors have been calculated and the bias in the sample calculations with the MCU code is approximated to Dkbias ¼ 0:00040. The largest effects on keff are caused by structures outside the reflectors (+0.00280) and the neutron detectors (-0.00111). Smaller effects are caused by impurities in the fuel particles and the CR6 control rod support structure with tubes all of which have not been modeled in the core. All other simplifications are considered to have a very limited effect. The impurities in the fuel particles were modeled in our model and of the core structures, only the CR6 support structure was not modeled. The overall difference from unity was approximated to 0.00140 based on the documentation and an assumption of how the reactivity effect is divided between the CR6 support structure and the experimental channels. Absolute values of the calculated multiplication factors were quite different depending on the material library used. With the three tested material libraries the difference between the highest and the lowest value in each case was nearly 2000 pcm. The cross sections of the graphite are believed to explain most of the differences between the libraries as graphite is the dominant material in the test facility. Different thermal scattering libraries JEF-2.2, ENDF/B-VI.8 and -VII, were also tested and each of them resulted in a slightly lower multiplication factor. The maximum difference was obtained with JEF-2.2, for which the difference was 0.00070. The offset of the averaged results from unity was then assumed to be caused by differences in the graphite absorption cross section. Sensitivity analysis that was performed with earlier model in Suikkanen et al. (2010) showed that difference of keff between JEFF-3.1.1, JEFF-3.1 and ENDF/B-VII were less than 0.001. This difference is very small if compared to difference between JEFF-3.1.1 and JENDL-4 in current results. In case of HTTR calculations similar difference between JENDL-4 and both JEFF-3.1 and ENDF/B-VII were reported by Shimakawa et al. (2010). Further, the difference in keff is explained to be caused by difference in graphite absorption cross section between the material libraries. MONK9A calculations of case 1 with BINGO material library are reported in appendix of benchmark documentation. The average result 1.0081 obtained with MONK9A for keff is very near to the JEFF-3.1.1 result of Serpent. Graphite in TENDL-2012 cases were modeled as C-12 and C-13 with natural isotopic fractions. To study the effect of graphite between TENDL-2012 and JENDL-4 the mixture of C-12 and C-13 was replaced by natural carbon from JENDL-4 library in case 3. After changing graphite the TENDL-2012 result was only 0.0005 larger than the JENDL-4 result. Graphite seems to cause the large differences between the results of different material libraries.

Another issue in the results was the growing trend of keff between different benchmark configurations in all the different calculation sets including the sample calculations provided with the benchmark documentation. The trend was very similar between the sample results and our results. In the sample calculations the pebble bed had been modeled with fuel elements in a regular lattice instead of a random packing as in our calculations. This trend is discussed briefly in the benchmark documentation where it is presented that the approximate account of the upper and lower part of the regular fuel element lattice causes a larger relative error in the case of a lower bed than with a higher bed. As the calculations with random packing show similar results when compared with the regular lattice, it is likely that pebble bed modeling is not the reason for the growing trend in the results. The original explanation might still hold if the local packing fraction variation were wrong similarly in both cases. This scenario is, however, considered unlikely. The pebble bed in MCU sample calculations has been generated by first creating densest possible packing of pebbles. Horizontal square lattice of pebbles is constructed and then these layers are piled on top of each other with making half pebble shift in both horizontal directions. After horizontal shifts the layer is moved upward until pebbles rest on top of the preceding layer. To obtain desired packing fraction the pitch of square lattice and the distance of layers are increased. To model individual fuel pebbles the coating layers of fuel particles are homogenized with fuel pebble graphite and positions for fuel particle kernels are generated on the fly during the neutron random walk (Ponomarev-Stepnoi et al., 2004). Various sensitivity analyses were made to find the source for the growing trend in the results. The previously mentioned material libraries and different thermal scattering libraries were tested and then tests were continued with and without the boron equivalent in graphite, impurities in the fuel particles, various Serpent related calculation options and even environmental conditions. The results of these analyses were consistent with the benchmark documentation and affected all cases equally enough, not to be the cause for the seen trend. The effect of the CR6 support structure was not calculated but some effect might be possible as the object is located in the core where it changes the packing of the pebbles in the bottom part of the pebble bed and additionally captures neutrons. 6. Conclusions The explicit random placement of fuel particles inside the pebbles and the use of DEM simulations to packing the pebbles add realism for the reactor physics calculations as they represent the actual geometric conditions and thus do not introduce additional modeling uncertainties, such as neutron streaming paths or intersected fuel particles and pebbles at domain boundaries. The use of DEM to construct pebble beds requires some additional effort but on the other hand similar effort might be spent on constructing equivalent regular lattice representation of the pebble bed. A method allowing exact pebble bed modeling available in Serpent was documented and demonstrated with the benchmark cases of ASTRA criticality experiments. Multiplication factors


V. Rintala et al. / Annals of Nuclear Energy 77 (2015) 223–230

increase with the height of the pebble bed in a similar way as in the sample results calculated with MCU Monte Carlo code. The absolute values are strongly dependent on the material library selected for the calculation. After a thorough search for errors we are confident that the HTR specific geometry types in Serpent are reliable and the deviation in the results is due to other factors. Serpent allows accurate Monte Carlo reactor physics modeling of HTR geometries with only a minor increase in calculation time compared to the conventional regular lattice approach. Serpent can also generate pebble-wise power distribution as output which is useful for thermal-hydraulics calculations. In further work, Serpent is planned to be used together with a thermal-hydraulic solver to obtain a coupled solution of reactor physics and thermalhydraulics. The DEM pebble packing simulation technique will be used to provide detailed data also for the thermal-hydraulic model. Acknowledgements This work was funded by the Academy of Finland under Grant 258145. We also acknowledge the work done at the RRC Kurchatov Institute for the ASTRA benchmark. References Brown, F.B., Martin, W.R., 2004. Stochastic geometry capability in MCNP5 for the analysis of particle fuel. Ann. Nucl. Energy 31, 2039–2047. Cundall, P.A., Strack, O.D.L., 1979. A discrete numerical model for granular assemblies. Géotechnique 29, 47–65. Di Renzo, A., Di Maio, F.P., 2004. Comparison of contact-force models for the simulation of collisions in DEM-based granular flow codes. Chem. Eng. Sci. 59, 525–541. Difilippo, F.C., 2003. Monte Carlo calculations of pebble bed benchmark configurations of the PROTEUS facility. Nucl. Sci. Eng. 143, 240–253. Dullien, F.A.L., 1979. Porous Media: Fluid Transport and Pore Structure. Academic Press, Inc., New York, NY. Garin, V., Glushkov, A., Glushkov, E., Gomin, E., Gurevich, M., Zimin, A., Kompaniets, G., Kukharkin, N., Lobyntsev, V., Nosov, V., Polyakov, D., Ponomarev-Stepnoi, N., Smirnov, O., Tel’kovskaya, O., Chunyaev, E., 2010. Evaluated benchmark experiments at critical assemblies simulating features of an HTHR at the ASTRA facility. Phys. At. Nucl. 73, 2271–2289.

Jodrey, W.S., Tory, E.M., 1985. Computer simulation of close random packing of equal spheres. Phys. Rev. A 32, 2347–2351. Kloss, C., Goniva, C., Hager, A., Amberger, S., 2012. Models, algorithms and validation for opensource DEM and CFD-DEM. Prog. Comput. Fluid Dynam. Int. J. 12, 140– 152. Leppänen, J., 2007. Randomly dispersed particle fuel model in the PSG Monte Carlo neutron transport code. In: M&C + SNA 2007, Monterey, CA. Leppänen, J., 2010. Performance of woodcock delta-tracking in lattice physics applications using the Serpent Monte Carlo reactor physics burnup calculation code. Ann. Nucl. Energy 37, 715–722. Melese, G., Katz, R., 1984. Thermal and Flow Design of Helium-cooled Reactors. American Nuclear Society, La Grange Park, Illinois. Onoda, G.Y., Liniger, E.G., 1990. Random loose packings of uniform spheres and the dilatancy onset. Phys. Rev. Lett. 64, 504–507. Plimpton, S., 1995. Fast parallel algorithms for short-range molecular dynamics. J. Comput. Phys. 117, 1–19. Ponomarev-Stepnoi, N.N., Bryzgalov, V.I., Glushkov, E.S., Gomin, E.A., Gurevich, M.I., Demin, V.E., Kompaniets, G.V., Lobyntsev, V.A., Nosov, V.I., Polyakov, D.N., Smirnov, O.N., Telk´ovskaya, O.V., 2004. Using the MCU computer program to analyze the results of critical experiments with HTGR fuel Pellets on the ASTRA testing stand. At. Energ. 97, 669–677. Ponomarev-Stepnoi, N.N., Glushkov, E.S., Kompaniets, G.V., Polyakov, D.N., 2008. Graphite annular core assemblies with spherical fuel elements containing coated UO2 fuel particles. In: International Handbook of Evaluated Reactor Physics Benchmark Experiments, March 2012 Edition. OECD NEA Nuclear Science Committee. NEA/NSC/DOC(2006)01. Rycroft, C.H., Grest, G.S., Landry, J.W., Bazant, M.Z., 2006. Analysis of granular flow in a pebble-bed nuclear reactor. Phys. Rev. E 74, 021306. Shimakawa, S., Goto, M., Nakagawa, S., Tachibana, Y., 2010. Impact of capture crosssection of carbon on nuclear design for HTGRs. In: Proceedings of the 5th International Conference on High Temperature Reactor Technology, HTR 2010, Prague, Czech Republic, October 18–20. Song, C., Wang, P., Makse, H.A., 2008. A phase diagram for jammed matter. Nature 453, 629–632. Suikkanen, H., Rintala, V., Kyrki-Rajamäki, R., 2010. An approach for detailed reactor physics modelling of randomly packed pebble beds. In: Proceedings of the 5th International Conference on High Temperature Reactor Technology, HTR 2010, Prague, Czech Republic, October 18–20. Suikkanen, H., Ritvanen, J., Jalali, P., Kyrki-Rajamäki, R., 2014. Discrete element modelling of pebble packing in pebble bed reactors. Nucl. Eng. Des. 273, 24–32. Visscher, W.M., Bolsterli, M., 1972. Random packing of equal and unequal spheres in two and three dimensions. Nature 239, 2727–2730. Wu, Z., Lin, D., Zhong, D., 2002. The design features of the HTR-10. Nucl. Eng. Des. 218, 25–32. Xiaowei, L., Xiaotian, L., Suyuan, Y., 2010. Nuclear graphite friction properties and the influence of friction properties on the pebble bed. Nucl. Eng. Des. 240, 2674– 2681.