This work is copyright 2010 by Joshua J. Cogliati. This work is licensed under the Creative Commons Attribution 3.0 Unported License. To view a copy of this license, visit http://creativecommons.org/licenses/by/3.0/ or send a letter to Creative Commons, 171 Second Street, Suite 300, San Francisco, California, 94105, USA. As specified in the license, verbatim copying and use is freely permitted, and properly attributed use is also allowed.

To the graduate faculty:

The members of the Committee appointed to examine the thesis of Joshua J. Cogliati find it satisfactory and recommend that it be accepted.

Thanks are due to many people who have provided information, comments and insight. I apologize in advance for anyone that I have left out. The original direction and technical ideas came from my adviser Abderrafi Ougouag as well as continued encouragement and discussion throughout its creation. Thanks go to Javier Ortensi for the encouragement and discussion as he figured out how use the earthquake data and I figured out how to generate it. At INL the following people assisted: Rob Bratton and Will Windes with the graphite literature review, Brian Boer with discussion and German translation, Hongbin Zhang with encouragement and Chinese translation and Suzette J. Payne for providing me with the Earthquake motion data. The following PBMR (South Africa) employees provided valuable help locating graphite literature and data: Frederik Reitsma, Pieter Goede and Alastair Ramlakan. The Juelich (Germany) people Peter Pohl, Johannes Fachinger and Werner von Lensa provided valuable assistance with understanding AVR. Thanks to Professor Mary Dunzik-Gougar for introducing me to many of these people, as well as encouragement and feedback on this PhD and participating as co-chair on the dissertation committee. Thanks to the other members of my committee, Dr. Michael Lineberry, Dr. Steve C. Chiu and Dr. Steve Shropshire, for providing valuable feedback on the dissertation. Thanks to Gannon Johnson for pointing out that length needed to be tallied separately from the length times force tally for the wear calculation (this allowed the vibration issue to be found). Thanks to Professor Jan Leen Kloosterman of the Delft University of Technology for providing me the PEBDAN program used for calculating Dancoff factors. The work was partially supported by the U.S. Department of Energy, Assistant Secretary for the office of Nuclear Energy, under DOE Idaho Operations Office Contract DEAC07-05ID14517. The financial support is gratefully acknowledged.

This dissertation contains work that was first published in the following conferences: Mathematics and Computation, Supercomputing, Reactor Physics and Nuclear and Biological Applications, Palais des Papes, Avignon, France, September 12-15, 2005; HTR2006: 3rd International Topical Meeting on High Temperature Reactor Technology October 1-4, 2006, Johannesburg, South Africa; Joint International Topical Meeting on Mathematics & Computation and Supercomputing in Nuclear Applications (M&C + SNA 2007) Monterey, California, April 15-19, 2007; Proceedings of the 4th International Topical Meeting on High Temperature Reactor Technology, HTR2008, September 28-October 1, 2008, Washington, DC USA and PHYSOR 2010 - Advances in Reactor Physics to Power the Nuclear Renaissance Pittsburgh, Pennsylvania, USA, May 9-14, 2010.

My mother helped me edit papers I wrote in my pre-college years, and my father taught me that math is useful in the real world. For this and all the other help in launching me into the world, I thank my parents.

Last, but certainly not least, thanks go to my wife, Elizabeth Cogliati for her encouragement and support. This support has provided me with the time needed to work on the completion of this dissertation. This goes above and beyond the call of duty since I started a PhD the same time we started a family. Thanks to my son for letting my skip my first Nuclear Engineering exam by being born, as well as asking all the really important questions. Thanks to my daughter for working right beside me on her dissertation on the little red toy laptop.

Contents

List of Figures

List of Tables

Abstract

Pebble bed reactors (PBR) have moving graphite fuel pebbles. This unique feature provides advantages, but also means that simulation of the reactor requires understanding the typical motion and location of the granular flow of pebbles.

This dissertation presents a method for simulation of motion of the pebbles in a PBR. A new mechanical motion simulator, PEBBLES, efficiently simulates the key elements of motion of the pebbles in a PBR. This model simulates gravitational force and contact forces including kinetic and true static friction. It’s used for a variety of tasks including simulation of the effect of earthquakes on a PBR, calculation of packing fractions, Dancoff factors, pebble wear and the pebble force on the walls. The simulator includes a new differential static friction model for the varied geometries of PBRs. A new static friction benchmark was devised via analytically solving the mechanics equations to determine the minimum pebble-to-pebble friction and pebble-to-surface friction for a five pebble pyramid. This pyramid check as well as a comparison to the Janssen formula was used to test the new static friction equations.

Because larger pebble bed simulations involve hundreds of thousands of pebbles and long periods of time, the PEBBLES code has been parallelized. PEBBLES runs on shared memory architectures and distributed memory architectures. For the shared memory architecture, the code uses a new O(n) lock-less parallel collision detection algorithm to determine which pebbles are likely to be in contact. The new collision detection algorithm improves on the traditional non-parallel O(n log(n)) collision detection algorithm. These features combine to form a fast parallel pebble motion simulation.

The PEBBLES code provides new capabilities for understanding and optimizing PBRs. The PEBBLES code has provided the pebble motion data required to calculate the motion of pebbles during a simulated earthquake. The PEBBLES code provides the ability to determine the contact forces and the lengths of motion in contact. This information combined with the proper wear coefficients can be used to determine the dust production from mechanical wear. These new capabilities enhance the understanding of PBRs, and the capabilities of the code will allow future improvements in understanding.

Chapter 1
Introduction

1.1 Pebble Bed Reactors Introduction

Pebble bed nuclear reactors are a unique reactor type that have been proposed and used experimentally. Pebble bed reactors were initially developed in Germany in the 1960s when the AVR demonstration reactor was created. The 10 megawatt HTR-10 reactor achieved first criticality in 2000 in China and future reactors are planned. In South Africa, Pebble Bed Modular Reactor Pty. Ltd. was designing a full scale pebble bed reactor to produce process heat or electricity.

Pebble bed nuclear reactors use graphite spheres (usually about 6 cm in diameter) for containing the fuel of the reactor. The graphite spheres encase smaller spheres of TRistructural-ISOtropic (TRISO) particle fuel. Unlike most reactors, the fuel is not placed in an orderly static arrangement. Instead, the graphite spheres are dropped into the top of the reactor, travel randomly down through the reactor core, and are removed from the bottom of the reactor. The pebbles are then possibly recirculated depending on the amount of burnup of the pebble and the reactor’s method of operation.

The first pebble bed reactor was conceived in 1950s in the West Germany using helium gas-cooling and spherical graphite fuel elements. Construction on the Arbeitsgemeinschaft Versuchsreaktor (AVR) 15 MWe reactor was started in 1959 at the KFA Research Centre Jülich. It started operation in 1967 and continued operation for 21 years until 1988. The reactor operated with an outlet temperature of 950$\circ$C. The AVR demonstrated the potential for the pebble bed reactor concept. Over the course of its operation, loss-of-coolant experiments were successfully performed.

The second pebble bed reactor was the Thorium High Temperature Reactor (THTR). This reactor was built in West Germany for an electric utility. It was a 300 MWe plant that achieved full power in September 1986. In October 1988, when the reactor was shutdown for maintenance, 35 bolt heads were found in the hot gas ducts leading to the steam generators. The determination was made that the plant could be restarted, but funding difficulties prevented this from occurring and the reactor was decommissioned (Goodjohn, 1991).

The third pebble bed reactor to be constructed and the only one that is currently operable is the 10 MWt High Temperature Reactor (HTR-10). This reactor is at the Tsinghua University in China. Construction was started in 1994 and reached first criticality in December 2000. This reactor is helium cooled and has an outlet temperature of 700$\circ$C (Wu et al., 2002; Xu and Zuo, 2002).

The use of high temperature helium cooled graphite moderated reactors with TRISO fuel particles have a number of advantages. A TRISO particle consists of spherical fuel kernel (such as uranium oxide) surrounded by four concentric layers: 1) a porous carbon buffer layer to accommodate fission product gases which limits pressure on the outer layers, 2) an interior pyrolytic carbon layer, 3) a layer of silicon carbide, and 4) an outer layer of pyrolytic carbon. The pyrolytic layers shrink and creep with irradiation, partially offsetting the pressure from the fission products in the interior as well as helping contain the fission gases. The silicon carbide acts as a containment mechanism for the metallic fission products.(Miller et al., 2002) These layers provide an in-core containment structure for the radioactive fuel and fission products.

The high temperature gas reactors have some advantages over conventional light water reactors. First, the higher outlet temperatures allow higher Carnot efficiency to be obtained¹ . Second, the higher temperatures can be used for process heat, which can reduce the use of methane. Third, the high temperature under which TRISO particles can operate allows for the exploitation of the negative temperature coefficient to safely shutdown the reactor without use of control rods.² Fourth, the higher temperature is above the annealing temperature for graphite, which safely removes Wigner energy³ . These are advantages of both prismatic and pebble bed high temperature reactors.(Gougar et al., 2004; Wu et al., 2002)

Pebble bed reactors, unlike most other reactors types, have moving fuel. This provides advantages but complicates modeling the reactors. A key advantage is that pebble bed reactors can be refueled online, that is, reactor shutdown is not needed for refueling. As a consequence, the reactors have low excess reactivity, as new pebbles can be added or excess pebbles removed to maintain the reactor at critical. The low excess reactivity removes the need for burnable poisons. A final advantage is that the moving fuel allows the pebble bed to be run with optimal moderation, where both increases and decreases in the fuel-to-moderator ratio cause reduction in reactivity. Ougouag et al. (2004) discuss the advantages of optimal moderation including improved fuel utilization. However, because the fuel is moving, many traditional methods of modeling nuclear reactors are inapplicable without a method for quantifying the motion. Hence, there is a need for development of methods usable for pebble bed reactor modeling.

1.2 Dissertation Introduction

This dissertation describes a computer code, PEBBLES, that is designed to provide a method of simulating the motion of the pebbles in a pebble bed reactor.

Chapter 4 provides the details of how the simulation works. Chapter 5 has a new static friction model developed for this dissertation.

Several checks have been made of the code. Figure 3.1 compares the PEBBLES simulation to experimentally determined radial packing fractions. Section 5.1 describes a new analytical benchmark that was used to test the static friction model in PEBBLES. Section 5.2 uses the Janssen model to test the static friction in a cylindrical vat.

Motivating all the above are the new applications, including Dancoff factors (8.1.1), calculating the angle of repose (8.1.2) and modeling an earthquake in section 8.2.

Chapter 2
Motivation

Most nuclear reactors have fixed fuel, including typical light water reactors. Some reactor designs, such as non-fixed fuel molten salt reactors, have fuel that is in fluid flow. Most designs for pebble bed reactors, however, have moving granular fuel. Since this fuel is neither fixed nor easily treatable as a fluid, predicting the behavior of the reactor requires the ability to understand the characteristics of the positions and motion of the pebbles. For example, predicting the probability of a neutron leaving one TRISO’s fueled region and entering another fueled region depends on the typical locations of the pebbles. A second example is predicting the effect of an earthquake on the reactivity of the pebble bed reactor. This requires knowing how the positions of the pebbles in the reactor change from the forces of the earthquake. Accurate prediction of the typical features of the flow and arrangement of the pebbles in the pebble bed reactor would be highly useful for their design and operation.

The challenge is to gain the ability to predict the pebble flow and pebble positions for start-up, steady state and transient pebble bed reactor operation.

The objective of the research presented in this dissertation is to provide this predicting ability. The approach used is to create a distinct element method computer simulation. The simulation determines the locations and velocities of all the pebbles in a pebble bed reactor and can calculate needed tallies from this data. Over the course of creating this simulation, various applications of the simulation were performed. These models allow the operation of the pebble bed reactor to be better understood.

Chapter 3
Previous work

Because the purpose of this dissertation is to produce a high fidelity simulation that can provide predictions of the pattern and flow of pebbles, previous efforts to simulate granular methods and packing were examined. A variety of simulations of the motion of discrete elements have been created for different purposes. Lu et al. (2001) applied a discrete element method (DEM) to determine the characteristics of packed beds used as fusion reactor blankets. Jullien et al. (1992) used a DEM to determine packing fractions for spheres using different non-motion methods. Soppe (1990) used a rain method to determine pore structures in different sized spheres. The rain method randomly chooses a horizontal position, and then lowers a sphere down until it reaches other existing spheres. This is then repeated to fill up the container. Freund et al. (2003) used a rain method for fluid flow in chemical processing.

The use of non-motion pebble packing methods provide an approximation of the positions of the pebble. Unfortunately, non-motion methods will tend to either under pack or over pack (sometimes both in the same model). For large pebble bed reactors, the approximately ten-meter height of the reactor core will result in different forces at the bottom than at the top. This will change the packing fractions between the top and the bottom, so without key physics, including static friction and the transmittal of force, non-motion physics models will not even be able to get correct positional information. Non-physics based modeling can not be used for predicting the effect of changes in static friction or pebble loading methods even if only the position data is required.

The initial PEBBLES code for calculation of pebble positions minimized the sum of the gravitational and Hookes’ law potential energies by adjusting pebble positions. However, that simulation was insufficient for determining flow and motion parameters and simulation of earthquake packing.

Additional references addressing full particle motion simulation were evaluated. Kohring (1995) created a 3-D discrete element method simulation to study diffusional mixing and provided detailed information on calculating the kinetic forces for the simulation. The author describes a simple method of calculating static friction. Haile (1997) discusses both how to simulate hard spheres and soft spheres using only potential energy. The soft sphere method in Haile proved useful for determining plausible pebble positions, but is insufficient for modeling the motion. Hard spheres are simulated by calculating the collision results from conservation laws. Soft spheres are simulated by allowing small overlaps, and then having a resulting force dependent on the overlap. Soft spheres are similar to what physically happens, in that the contact area distorts, allowing distant points to approach closer than would be possible if the spheres were truly infinitely hard and only touched at one infinitesimal point. Hard spheres are impractical for a pebble bed due to the frequent and continuous contact between spheres so soft spheres are used instead. The dissertation by Ristow (1998) describes multiple methods for simulation of granular materials. On Ristow’s list of methods was a model similar to that used as the kernel of the work supporting this dissertation. Ristow’s dissertation mentioned static friction and provided useful references that will be discussed in Section 3.2.

To determine particle flows, Wait (2001) developed a discrete element method that included only dynamic friction. Concurrently with this dissertation research, Rycroft et al. (2006b) used a discrete element method, created for other purposes, to simulate the flow of pebbles through a pebble bed reactor.

Multiple other discrete element codes have been created and PEBBLES is similar to several of the full motion models. For most of the applications discussed in this dissertation, only a model that simulates the physics with high fidelity is useful. The PEBBLES dynamic friction model is similar to the model used by Wait or Rycroft, but the static friction model incorporates some new improvements that will be discussed later.

In addition to simulation by computer, other methods of determining the properties of granular fluids have been used. Bedenig et al. (1968) used a scale model to experimentally determine residence spectra (the amount of time that pebbles from a given group take to pass through a reactor) for different exit cone angles. Kadak and Bazant (2004) used scale models and small spheres to estimate the flow of pebbles through a full scale pebble bed reactor. These researchers also examined the mixing that would occur between different radial zones as the pebbles traveled downward. Bernal et al. (1960) carefully lowered steel spheres into cylinders and shook the cylinders to determine both loose and dense packing fractions. The packing fraction and boundary density fluctuations were experimentally measured by Benenati and Brosilow (1962). The Benenati and Brosilow data have been used to verify that the PEBBLES code was producing correct boundary density fluctuations (See Figure 3.1). Many experiments were performed in the designing and operating of the AVR reactor to determine relevant properties such as residence times and optimal chute parameters (Bäumer et al., 1990). These experiments provide data for testing the implementation of any computational model of pebble flow.

The PEBBLES simulation uses elements from a number of sources and uses standard classical mechanics for calculating the motion of the pebbles based on the forces calculated. The features in PEBBLES have been chosen to implement the necessary fidelity required while allowing run times small enough to accommodate hundreds of thousands of pebbles. The next sections will discuss handling static friction.

3.1 Static Friction Overview

Static friction is an important effect in the movement of pebbles and their locations in pebble bed reactors. This section briefly reviews static friction and its effects in pebble bed reactors. Static friction is a force between two contacting bodies that counteracts relative motion between them when they are moving sufficiently slowly(Marion and Thornton, 2004). Macroscopically, the maximum magnitude of the force is proportional to the normal force with the following equation:

where $\mu$ is the coefficient of static friction, $F_s$ is the static friction force and $F_{\perp }$ is the normal (load) force.

Static friction results in several effects on granular materials. Without static friction, the angle of the slope of a pile of a material (angle of repose) would be zero(Duran, 1999). Static friction also allows ‘bridges’ or arches to be formed near the outlet chute. If the outlet chute is too small, the bridging will be stable enough to clog the chute. Static friction will also transfer force from the pebbles to the walls. This will result in lower pressure on the walls than would occur without static friction(Sperl, 2006; Walker, 1966).

For an elastic sphere, static friction’s counteracting force is the result of elastic displacement of the contact point. Without static friction, the contact point would slide as a result of relative motion at the surface. With static friction, the spheres will experience local shear that distorts their shape so that the contact point remains constant. This change will be called stuck-slip, and continues until the counteracting force exceeds $\mu F_{\perp }$. When the counteracting force exceeds that value, the contact point changes and slide occurs. The mechanics of this force with elastic spheres were investigated by Mindlin and Deresiewicz (1953). Their work created exact formulas for the force as a function of the past relative motion and force.

3.2 Simulation of Mechanics of Granular Material

Many simulations of granular materials incorporating static friction have been devised. Cundall and Strack (1979) developed an early distinct element simulation of granular materials that incorporated a computationally efficient static friction approximation. Their method involved integration of the relative velocity at the contact point and using the sum as a proxy for the current static friction force. Since their method was used for simulation of 2-D circles, adaptation was required for 3-D granular materials. One key aspect of adaptation is determining how the stuck-slip direction changes as a result of contacting objects’ changing orientation.

Vu-Quoc and Zhang (1999) and Vu-Quoc et al. (2000) developed a 3-D discrete-element method for granular flows. This model was used for simulation of particle flow in chutes. They used a simplification of the Mindlin and Deresiewicz model for calculating the stuck-slip magnitude, and project the stuck-slip onto the tangent plane each time-step to rotate the stuck-slip force direction. This correctly rotates the stuck-slip, but requires that the rotation of the stuck-slip be done as a separate step since it not written in a differential form.

Silbert et al. (2001) and Landry et al. (2003) describe a 3-D differential version of the Cundall and Strack method. The literature states that particle wall interactions are done identically. The amount of computation of the model is less than the Vu-Quoc, Zhang and Walton model. This model was used for modeling pebble bed flow(Rycroft et al., 2006a,b). This model however, does not specify how to apply their differential version to modeling curved walls.

Chapter 4
Mechanics Model

The PEBBLES simulation calculates the forces on each individual pebble. These forces are then used to calculate the subsequent motion and position of the pebbles.

4.1 Overview of Model

The PEBBLES simulation tracks each individual pebble’s velocity, position, angular velocity and static friction loadings. The following classical mechanics differential equations are used for calculating the time derivatives of those variables:

where $F_{ij}$ is the force from pebble $j$ on pebble $i$, $F_{ci}$ is the force of the container on pebble $i$, $g$ is the gravitational acceleration constant, $m_i$ is the mass of pebble $i$, $v_i$ is the velocity of pebble $i$, $p_i$ is the position vector for pebble $i$, ${\omega }_i$ is the angular velocity of pebble i, $F_{\parallel ij}$ is the tangential force between pebbles $i$ and $j$, $F_{\perp ij}$ is the perpendicular force between pebbles $i$ and $j$, $r_i$ is the radius of pebble $i$, $I_i$ is the moment of inertia for pebble i, $F_{\parallel ci}$ is the tangential force of the container on pebble i, ${\widehat{n}}_{ci}$ is the unit vector normal to the container wall on pebble i, ${\widehat{n}}_{ij}$ is the unit vector pointing from the position of pebble $i$ to that of pebble $j$, $s_{ij}$ is the current static friction loading between pebbles $i$ and $j$, and $S$ is the function to compute the change in the static friction loading. The static friction model contributes to the $F_{\parallel ij}$ term which is also part of the $F_{ij}$ term. Figure 4.1 illustrates the principal vectors with pebble $i$ going in the $v_i$ direction and rotating around the ${\omega }_i$ axis, and pebble $j$ going in the $v_j$ direction and rotating around the ${\omega }_j$ axis.

The mass and moment of inertia are calculated assuming spherical symmetry with the equations:

where $r_c$ is the radius of inner (fueled) zone of the pebble, $r_o$ is the radius of whole pebble, ${\rho }_c$ is the average density of center fueled region and ${\rho }_o$ is the average density of outer non-fueled region.

The dynamic (or kinetic) friction model is based on the model described by Wait (2001). Wait’s and PEBBLES model calculate the dynamic friction using a combination of the relative velocities and pressure between the pebbles, as shown in Equations (4.7) and (4.8):

where $C_{\parallel }$ is the tangential dash-pot constant, $C_{\perp }$ is the normal dash-pot constant, $F_{\perp ij}$ is the normal force between pebbles i and j, $F_{d\parallel ij}$ is the tangential dynamic friction force between pebbles i and j, $h$ is the normal Hooke’s law constant, $l_{ij}$ is the overlap between pebbles i and j, $v_{\parallel ij}$ is the component of the velocity between two pebbles perpendicular to the line joining their centers, $v_{\perp ij}$ is the component of the velocity between two pebbles parallel to the line joining their centers, $v_{ij}$ is the relative velocity between pebbles i and j and $\mu$ is the kinetic friction coefficient. Equations (4.9-4.12) relate supplemental variables to the primary variables:

The friction force is then bounded by the friction coefficient and the normal force, to prevent it from being too great:

where $F_{s\parallel ij}$ is the static friction force between pebbles $i$ and $j$, $F_{d\parallel ij}$ is the kinetic friction force between pebbles $i$ and $j$, $h_s$ is the coefficient for force from slip, $s_{ij}$ is the slip distance perpendicular to the normal force between pebbles $i$ and $j$, $v_{max}$ is the maximum velocity under which static friction is allowed to operate, and $\mu$ is the static friction coefficient when the velocity is less than $v_{max}$ and the kinetic friction coefficient when the velocity is greater. These equations fully enforces the first requirement of a static friction method, $\left|F_s\right|\le \mu \left|F_{\perp }\right|$.

4.2 Integration

When all the position, linear velocity, angular velocity and slips are combined into a vector $y$, the whole computation can be written as a differential formulation in the form:

This can be solved by a variety of methods with the simplest being Euler’s method:

In addition, both the Runge-Kutta method and the Adams-Moulton method can be used for solving this equation. These methods improve the accuracy of the simulation. However, they do not improve the wall-clock time at the lowest stable simulation, since the additional time required for computation negates the advantage of being able to use somewhat longer time-steps. In addition, when running on a cluster, more data needs to be transferred since the methods allow non-contacting pebbles to affect each other in a single ‘time-step calculation’.

4.3 Geometry Modeling

For any geometry interaction, two things need to be calculated, the overlap distance $l$ (or, technically, the mutual approach of distant points) and the normal to the surface $\widehat{n}$. The input is the radius of the pebble $r$ and the position of the pebble, $p$ with components $p_x$, $p_y$, and $p_z$

For the floor contact this is:

For cylinder contact on the inside of a cylinder this is:

For cylinder contact on the outside of a cylinder this is:

For contact on the inside of a cone defined by the $radius=mz+b$:

These equations are derived from minimizing the distance between the contact point $c$ and the pebble position $p$.

For contact on a plane defined by $ax+by+cz+d=0$ where the equation has been normalized so that $a^2+b^2+c^2=1$, the following is used:

Combinatorial geometry operations can be done. Intersections and unions of multiple geometry types are done by calculating the overlaps and normals for all the geometry objects in the intersection or union. For an intersection, where there is overlap on all the geometry objects, then the smallest overlap and associated normal are kept, which may be no overlap. For a union, the largest overlap and its associated normal are kept.

For testing that a geometry is correct, a simple check is to fill up the geometry with pebbles using one of the methods described in Section 4.4, and then make sure that linear and angular energy dissipate. Many geometry errors will show up by artificially creating extra linear momentum. For example, if a plane is only defined at the top, but it is possible for pebbles to leak deep into the bottom of the plane, they will go from having no overlap to a very high overlap, which will give the pebble a large force. This results in extra energy being added each time a pebble encounters the poorly defined plane, which will show up in energy tallies.

4.4 Packing Methods

The pebbles are packed using three main methods. The simplest creates a very loose packing with an approximately 0.15 packing fraction by randomly choosing locations, and removing the overlapping ones. These pebbles then allowed to fall down to compact to a realistic packing fraction.

The second is the PRIMe method developed by Kloosterman and Ougouag (2005). In this method large numbers of random positions (on the order of 100,000 more than will fit) are generated. The random positions are sorted by height, and starting at the bottom, the ones that fit are kept. Figure 4.2 illustrates this process. This generates packing fractions of approximately 0.55. Then they are allowed to fall to compact. This compaction takes less time than starting with a 0.15 packing fraction.

The last method is to automatically generates virtual chutes above the bed where the actual inlet chutes are, and then loads the pebbles into the chutes. Figure 4.3 shows this in progress. This allows locations that have piles where the inlet chutes are, but can be done quicker than a recirculation. The other two methods generate flat surfaces at the top, which is unrealistic, since the surface of a recirculated bed will have cones under each inlet chute.

4.5 Typical Parameters

The typical parameters used with the PEBBLES code are described in Table 4.1. Alternative numbers will be described when used.

Chapter 5
A New Static Friction Model

The static friction model in PEBBLES is used to calculate the force and magnitude of the static friction force. Other models have been created before to calculate static friction, but the PEBBLES model provides the combination of being a differential model (as opposed to one where the force is rotated as a separate step) and being able to handle the type of geometries that exist in pebble bed reactors.

The static friction model has two key requirements. First, the force from stuck-slip must be updated based on relative motion of the pebbles. Second, the current direction of the force must be calculated since the pebbles can rotate in space.

5.0.1 Use of Parallel Velocity for Slip Updating

For elastic spheres, the true method of updating the stuck-slip force is to use the method of Mindlin and Deresiewicz (1953). This method requires computationally and memory intensive calculations to track the forces. Instead, a simpler method is used to approximate the force. This method, described by Cundall and Strack (1979) uses the integration of the parallel relative velocity as the displacement. The essential idea is that the farther the pebbles have stuck-slipped at the contact point, the greater the counteracting static friction force needs to be. This is what happens under more accurate models such as Mindlin and Deresiewicz. There are two approximations imposed by this assumption. First, the amount the force changes is independent of the normal force. Second, the true hysteretic effects that are dependent on details of the loading history are ignored. For simulations where the exact dynamics of static friction are important, these could potentially be serious errors. However, since static friction only occurs when the relative speed is low, the dynamics of the simulation usually are unimportant. Thus, for most circumstances, the following approximation can be used for the rate of change of the magnitude and non-rotational change of the stuck-slip:

5.0.2 Adjustment of Oversize Slips

The slips can build up to unrealistically large amounts, that is, greater than $\mu \left|F_{\perp }\right|$; equation 5.1 places no limit on the maximum size of the slip. The excessive slip is solved at two different locations. First, when the frictions are added together to determine the total friction they are limited by $\mu \left|F_{\perp }\right|$ in equation (4.14). This by itself is insufficient, because the slip is storing potential energy that appears anytime the normal force increases. This manifests itself by causing vibration of the pebbles to continue for long periods of time. Two methods for fixing the hidden slip problem are available in PEBBLES. The simplest drops any slip that exceeds the static friction threshold (or an input parameter value somewhat above the static friction threshold so small vibrations do not cause the slip to disappear).

The second method used in PEBBLES is to decrease the slip that is over a threshold value. If the slip is too great, a derivative that is the opposite as the current slip is added as an additional term in the slip time derivative. This occurs in the following additional term:

In this $R\left(x\right)$ is the ramp function (which is $x$ if $x>0$ and 0 otherwise), $s_{sd}$ is a constant to select how much the slip is allowed to exceed the static friction threshold (usually 1.1 in PEBBLES). This derivative adder is used in most PEBBLES runs since it does allow vibrational energy to decrease, yet does not cause the pyramid benchmark to fail like complete removal of too great slips does.

When using non-Euler integration methods, the change in slip is calculated multiple times. Each time it is calculated, it might be set to be zeroed. In the PEBBLES code, if any of the added up slips for a given contact were set to be zeroed, the final slip is zeroed. This is not an ideal method, but it works well enough.

5.0.3 Rotation of Stuck-Slip

The static friction force must also be rotated so that it is in the plane of contact between the two pebbles. When there is a difference between the pebbles’ center velocities, which changes in the relative pebble center location, change in the direction in the stuck-slip occurs. That is:

First, let $n_{ijn}=p_i-p_j$ and $d{n}_{ijn}=v_i-v_j$. The cross product $-d{n}_{ijn}\times n_{ijn}$ is perpendicular to both $n$ and $dn$ and signed to create the axis around which $s$ is rotated in a right-handed direction. Then, using the cross product of the axis and $s$, $-\left(d{n}_{ij}\times n_{ijn}\right)\times s_{ijn}$ gives the correct direction that $s$ should be increased.

Next, determine the factors required to make the differential the proper length. By cross product laws,

where $𝜃$ is the angle between $n_{ijn}$ and $d{n}_{ij}$ and $\varphi$ is the angle between $d{n}_{ij}\times n_{ijn}$ and $s_{ijn}$.

The relevant vectors are shown in figure 5.1.

The goal is to rotate $s$ by angle $\alpha \prime$ which is the ‘projection’ into the proper plane of the angle $\alpha$ that $n$ rotates by. Since the direction has been determined, for simplicity the figure leaves the indexes off, and concentrates on determining the lengths. In Figure 5.1, $s$ is the old slip vector, $s\prime$ is the new slip vector, $n$ is the vector pointing from one pebble to another. The vector $dn\Delta t$ is added to $n$ to get the new $n\prime$, $n+dn\Delta t$. The initial condition is that $s$ and $n$ are perpendicular. The final conditions are that $s\prime$ and $n\prime$ are perpendicular, and that $s$ and $s\prime$ are the same length and that $s\prime$ is the closest vector to $s$ as it can be while satisfying the other conditions. There is no requirement that $s$ or $s\prime$ are coplanar with $dn\Delta t$ (otherwise $\alpha \prime$ would be equal to $\alpha$). From the law of sines we have:

so

In Figure 5.2 the projection to the correct plane occurs. First by using $\varphi$ the length of $s$ is projected to the plane. Note that $\varphi$ is the angle both to $s$ and to $s\prime$. So, the length of the line on the $dn\times n$ plane is $\left|s\right|sin\varphi$, and the length of the straight line at the end of the triangle is $\left|s\right|sin\varphi sin\alpha$ (note that the chord length is $\left|s\right|\left(sin\varphi \right)\alpha$, but as $\Delta t$ approaches 0 the other can be used). From these calculations, the length of the $ds\Delta t$ can be calculated with the following formula:

Since $|-\left(d{n}_{ij}\times n_{ijn}\right)\times s_{ijn}|=\left|d{n}_{ij}\right|\left|n_{ijn}\right|\left|s_{ijn}\right|sin𝜃sin\varphi$ the formula for the rotation is:

As a differential equation this is:

By the vector property $a\times \left(b\times c\right)=b\left(a\cdot c\right)-c\left(a\cdot b\right)$ and since $\left(p_i-p_j\right)\cdot s_{ij}=0$, this can be rewritten as the version in Silbert et al. (2001):

5.0.4 Differential Equation for Surface Slip Rotating

It might seem that the wall interaction could be modeled the same way as the pebble-to-pebble interaction. For sufficiently simple wall geometries this may be possible, but actual pebble bed reactor geometries are more complicated, and violate some of the assumptions that underpin the derivation. For a flat surface, there is no rotation, so the formula can be entirely dropped. For a spherical surface, it would be possible to measure the curvature at pebble to surface contact point in the direction of relative velocity to the surface. This curvature could then be used as an effective radius in the pebble-to-pebble formulas.

The pebble reactor walls have additional features that violate assumptions made for the derivation. For surfaces such as a cone, the curvature is not in general constant, because the path can follow elliptical curves. As well, the curvature has discontinuities where different parts of the reactor join together (for example, the transition from the outlet cone to the outlet chute). At these transitions, the assumption that the slip stays parallel to the surface fails because the slip is parallel to the old surface, but the new surface has a different normal.

Because of the complications with using the pebble to pebble interaction, PEBBLES uses an approximation of the “rotation delta.” This is similar to the Vu-Quoc and Zhang (1999) method of adjusting the slip so that it is parallel to the surface each time. Each time when the slip is used, a temporary version of the slip that is properly aligned to the surface is computed and used for calculating the force. As well, a rotation to move the slip more parallel to the surface is also computed.

The rotation is computed as follows. Let the normal direction of the wall at the point of contact of the pebble be $n$, and the old stuck-slip be $s$. Let $a$ be the angle between $n$ and $s$. $n\times s$ is perpendicular to both $n$ and $s$ and so this cross product is the axis that needs to be rotated around. Then $\left(n\times s\right)\times s$ is perpendicular to this vector, so it is either pointing directly towards $n$ if $a$ is acute or directly away from $n$ if $a$ is obtuse. To obtain the correct direction, this vector is multiplied by the scalar $s\cdot n$ which has the correct sign from $cosa$. The magnitude of $\left(s\cdot n\right)\left[\left(n\times s\right)\times s\right]$ needs to be determined for reasonableness. We define the angle $b$, which is between $\left(n\times s\right)$ and $s$. By these definitions the magnitude is $\left(\left|s\right|\left|n\right|cosa\right)\left[\left(\left|n\right|\left|s\right|sina\right)\left|s\right|sinb\right]$. $b$ is a right angle since $n\times s$ is perpendicular to $s$, so $sinb=1$. Collecting terms gives the magnitude as $\left|s{|}^3\right|n{|}^{2}cosasina$ which is divided by $|n\times s|\left|n\right|\left|s\right|$ to give the full term the magnitude $\left|s\right|cosa$. This is the length of the vector that goes from $s$ to the plane perpendicular to $n$. This produces equation 5.11, which can be used to ensure that the wall stuck-slip vector rotates towards the correct direction.

5.1 Testing of Static Friction Model With Pyramid Test

Static friction is an important physical feature in the implementation of mechanical models of pebbles motion in a pebble bed, and checking its correctness is important. A pyramid static friction test model was devised as a simple tool for verifying the implementation of a static friction model within the code. The main advantages of the pyramid test are that the model test is realistic and that it can be modeled analytically, providing an exact basis for the comparison. The test benchmark consists of a pyramid of five spheres on a flat surface. This configuration is used because the forces acting on each pebble can be calculated simply and the physical behavior of a model with only kinetic friction is fully predictable on physical and mathematical grounds: with only kinetic friction and no static friction, the pyramid will quickly flatten. Even insufficient static friction will result in the same outcome. The four bottom spheres are arranged as closely as possible in a square, and the fifth sphere is placed on top of them as shown in Fig. 5.4.

The lines connecting the centers of the spheres form a pyramid with sides $2R$, as shown in Fig. 5.5, where $R$ is the radius of the spheres. The length of $a$ in the figure is $\frac{2R}{\sqrt{2}}$, and because $b$ is part of a right triangle, ${\left(2R\right)}^2-{\left(\frac{2R}{\sqrt{2}}\right)}^2=b^2=4R^2-\frac{4R^2}{2}=2R^2$, so $b$ has the same length as $a$, and thus the elevation angle for all vertexes of the pyramid are $45^{\circ }$ from horizontal.

Taking for origin of the coordinates system the projection of the pyramid summit onto the ground, the locations (coordinates) of the sphere centers are given in Table 5.1.

5.1.1 Derivation of Minimum Static Frictions

The initial forces on the base sphere are the force of gravity $mg$, and the normal forces $Tn$ and $Fn$ as shown in Fig. 5.6. This causes initial stuck-slip which will cause $Fs$ to develop to counter the slip, and $Ts$ to counter the rotation of the base sphere relative to the top sphere. The top sphere will have no rotation because the forces from the four spheres will be symmetric and counteract each other.

The forces on the base sphere are:

Note that $Fn$ is larger than $Tn$ since $Tn$ is only a portion of the $mg$ force since the top sphere transmits (and splits) its force onto all four base spheres.

There are three requirements for a base sphere to be non-accelerated.

If a base sphere is not rotating than there is no torque, so:

The resultant of all forces must also be zero in the x and the y direction (vector notation dropped since they are in one dimension and therefore scalars) as follows:

Since the angle of contact between a base sphere and the top sphere is $45^{\circ }$, the following two equations hold (where $Ts$ is the magnitude of $Ts$ and $Tn$ is the magnitude of $Tn$):

This changes equations 5.13 and 5.14 into:

Combining equation 5.12 and 5.17 provides:

Which gives the relation:

By the static friction Equation 3.1,

Combining equations 5.20 and 5.21 and simplifying gives the requirement that

For use with testing, the static friction program can be tested twice with a sphere-to-sphere friction coefficient slightly above 0.41421... and one slightly below 0.41421.... In the first case the pyramid should be stable, and in the second case the top ball should fall to the floor.

Since ¼ of the weight of the top pebble is on one of the base pebbles, the following holds:

Combining 5.18 and 5.23 provides the following equation:

Equations 5.17 and 5.24 can be added to produce

Using 5.12 and 5.24 and solving for Fs gives the following value for Fs:

By the static friction Equation 3.1:

Substituting the values for $Fs$ and $Fn$ gives:

Simplifying provides the following relation for the surface-to-sphere static friction requirement:

This can be used similarly to the other static friction requirement by setting the value slightly above 0.08284... and slightly below 0.08284... and making sure that it is stable with the higher value and not stable with the lower value.

This test was inspired by an observation of lead cannon balls stacked into a pyramid. I tried to stack used glass marbles into a five ball pyramid and it was not stable. Since lead has a static friction coefficient around 0.9 and used glass has a much lower static friction, the physics of pyramid stability was further investigated, resulting in this benchmark test of static friction modeling.

5.1.2 Use of Benchmark

The benchmark test of two critical static friction coefficients is defined by the following equations. If both static friction coefficients are above the critical values, the spheres will form a stable pyramid. If either or both values are below the critical values the pyramid will collapse.

To set up the test cases, the sphere locations from Table 5.1 should be used as the initial locations of the sphere. Then, static friction coefficients for the sphere-to-sphere contact and the sphere-to-surface contact are chosen. The code is then run until either the center sphere falls to the surface, or the pyramid obtains a stable state. There are three test cases that are run to test the model.

1. ${\mu }_{surface}={\mu }_{criticalsurface}+𝜖$ and ${\mu }_{sphere}={\mu }_{criticalsphere}+𝜖$ which should result in a stable pyramid.
2. ${\mu }_{surface}={\mu }_{criticalsurface}-𝜖$ and ${\mu }_{sphere}={\mu }_{criticalsphere}+𝜖$ which should result in a fall.
3. ${\mu }_{surface}={\mu }_{criticalsurface}+𝜖$ and ${\mu }_{sphere}={\mu }_{criticalsphere}-𝜖$ which should result in a fall.

For soft sphere models, there are fundamental limits to how close the model’s results can be to the critical coefficient. Essentially, as the critical coefficients are approached, the assumptions become less valid. For example, with soft (elastic) spheres, the force from the center sphere will distort the contact angle, so the actual critical value will be different. For the PEBBLES code, an $𝜖$ value of 0.001 is used.

5.2 Testing of the Static Friction Model by Comparison with Janssen’s Formula

The pyramid static friction test is used as a simple test of the static friction model. Another test compares the static friction model against the Janssen formula’s behavior (Sperl, 2006). This formula specifies the expected wall pressure as a function of depth. This formula only applies when the static friction is fully loaded, that is when $F_s|=\mu |F_{\perp }|$. This condition is generally not satisfied until some recirculation has occurred. Figure 5.7 shows the normal force and the static friction force from a pebble to the wall. With the PEBBLES code, this is only satisfied after recirculation with lower values of the static friction coefficient $\mu$.

The equation used to calculate the pressure on the region $R$ from the normal force in PEBBLES is:

where $p$ is the pressure, $R_h$ is the height of the region, and $r$ is the radius of the cylinder.

The equation for calculating the Janssen formula pressure is

where ${\mu }_{pp}$ is the pebble to pebble static friction coefficient, ${\mu }_{wall}$ is the pebble to wall, $f$ is the packing fraction, $\rho$ is the density, $g$ is the gravitational acceleration, and $z$ is the depth that the pressure is being calculated. For the Figures 5.8 and 5.9, a packing fraction of 0.61 is used and a density of 1760 kg/m${}^3$ are used. There are 20,000 pebbles packed into a 0.5 meter radius cylinder, and 1,000 are recirculated before the pressure measurement is done.

Figure 5.8 compares the Janssen model with the PEBBLES simulation for static friction values of 0.05 and 0.15. For this case, the Janssen formula and the simulated pressures match closely. Figure 5.9 compares these again. In this case, the 0.25 $\mu$ values only approximately match, and the 0.9 static friction pressure values do not match at all. The static friction slip vectors were examined, and they are not perfectly vertical, and they are not fully loaded. This results in the static friction force being less than the maximum possible, and thus the pressure is higher since less of the force is removed by the walls.

Chapter 6
Unphysical Approximations

Modeling the full physical effects that occur in a pebble bed reactor mechanics is not computationally possible with current computer resources. In fact, even modeling all the intermolecular forces that occur between two pebbles at sufficient levels to reproduce all macroscopic behavior is probably computationally intractable at the present time. This is partially caused by the complexity of effects such as inter-grain boundaries and small quantities of impurities that affect the physics and different levels between the atomic effects and the macroscopic world. Instead, all attempts at modeling the behavior of pebble bed reactor mechanics have relied on approximation to make the task computationally practical. The PEBBLES simulation has as high or higher fidelity than past efforts, but it does use multiple unphysical approximations. This chapter will discuss the approximations so that future simulation work can be improved, and an understanding of what limitations exist when applying PEBBLES to different problems.

In different regions of the reactor, the radioactivity and the fission will heat the pebbles differently, and the flow of the coolant helium will distribute this heat around the reactor. This will change the temperature of different parts of the reactor. Since the temperature will be different, the parameters driving the mechanics of the pebbles will be different as well. This includes parameters such as the static friction coefficients and the size of the pebbles which will change through thermal expansion. As well, parameters such as static friction can also vary depending on the gas in which they currently are in and in which they were, since some of the gas tends to remain in and on the carbon surface. Graphite dust produced by wear may also affect static friction in downstream portions of the reactor.

The pebbles in a pebble bed reactor have helium gas flowing around and past them. PEBBLES and all other pebble bed simulations ignore effects of this on pebble movement. However, the gas will cause both additional friction when the pebbles are dropping through the reactor, and the motion of gas will cause additional forces on pebbles.

Pebble bed mechanics simulations use soft spheres. Physically, there will be deflection of spheres under pressure (even the pressure of just one sphere on the floor), but the true compression is much smaller than what is actually modeled. In PEBBLES, the forces are chosen to keep the compression distance at a millimeter or below. Another effect related to the physics simulation is that force is transmitted via contact. This means the force from one end of the reactor is transmitted at a speed related to the time-step used for the simulation, instead of the speed of sound.

Since simulating billions of time-steps is time consuming, two approximations are made. First, instead of simulating the physical time that pebble bed reactors have between pebble additions (on the order of 2-5 minutes), new pebbles are added at a rate between a quarter second and two seconds. This may result in somewhat unphysical simulations since some vibration that would have dampened out with a longer time between pebble additions still exists when the next pebble impacts the bed. Second, since full recirculation of all the pebbles is computationally costly, for some simulations, only a partial recirculation or no recirculation is done.

The physics models do not take into account several physical phenomena. The physics do not handle pure spin effects, such as when two pebbles are contacting and are spinning with an axis around the contact point. This should result in forces on the pebbles, but the physics model does not handle this effect since the contact velocity is calculated as zero. In addition, when the pebble is rolling so that the contact velocity is zero because the pebble’s turning axis is parallel to the surface and at the same rate as the pebble is moving along the surface, there should be rolling friction, but this effect is not modeled either. As well, the equations used assume that the pebbles are spherically symmetric, but defects in manufacturing and slight asymmetries in the TRISO particle distribution mean that there will be small deviations from being truly spherically symmetric.

The physics model does not match classical Hertzian or Mindlin and Deresiewicz elastic sphere behavior. The static friction model is a simplification and does not capture all the hysteretic effects of true static friction. Effectively, this means that $h_s$, the coefficient used to calculate the force from slip, is not a constant. In order to fully discuss this, some features of these models will be discussed in the following paragraphs.

Since closed-form expressions exist for elastic contact between spheres, they will be used, instead of a more general case which lacks closed-form expressions. Spheres are not a perfect representation of the effect of contact between shapes such as a cone and a sphere, but should give an approximation of the size of the effect of curvature.

The amount of contact area and displacement of distant points for two spheres or one sphere and one spherical hole (that is negative curvature) for elastic spheres can be calculated via Hertzian theory(Johnson, 1985). For two spherical surfaces the following variables are defined:

and

with $R_i$ the $i$th’s sphere’s radius, $E_i$ the Young’s modulus, ${\nu }_i$ the Poisson’s ration of the material. For a concave sphere, the radius will be negative. Then, via Hertzian theory, the contact circle radius will be:

where $N$ is the normal force. The mutual approach of distant points is given by:

Notice that the above differs compared to the Hooke’s Law formulation that PEBBLES uses. The maximum pressure will be:

So as a function of the radii $R_1$ and $R_2$, the circle radius of the contact will be:

If $m$ is used for the multiple of negative curvature sphere of the radius of the other, then the equation becomes:

which can be rearranged to:

From this equation, as $m$ increases, it has less effect on contact area, so if $R_2$ is much greater than $R_1$, the contact area will tend to be dominated by $R_1$ rather than $R_2$. For example, typical radii in PEBBLES might be 18 cm outlet chute and a 3 cm pebble, which would put $m$ at 6, so the effect on contact area radius would be about 33% difference compared to pebble to pebble contact area radius, or 6% compared to a flat surface.¹

To some extent, the actual contact area is irrelevant for calculating the maximum static friction force as long as some conditions are met. Both surfaces need to be of a uniform material. The basic macroscopic description $\left|F_S\right|<=\mu \left|N\right|$ needs to hold, so changing the area changes the pressure $P=N∕a$, but not the maximum static friction force. If the smaller area causes the pressure to increase enough to cause plastic rather than elastic contact, then through that mechanism, the contact area would cause actual differences in experimental values. PEBBLES also does not calculate plastic contact effects.

The contact area causes an effect through another mechanism. The tangential compliance in the case of constant normal and increasing tangential force, that is the slope of the curve relating displacement to tangential force, is given in Mindlin and Deresiewicz as:

Since the contact area radius, $a$, is a function of curvature, the slope of the tangential compliance will be as well, which is another effect that PEBBLES’ constant $h_s$ does not capture.

In summary for the static friction using a constant coefficient for $h_s$ yields two different approximations. First, using the same constants for wall contact when there is different curvatures is an approximation that will give somewhat inconsistent results. Since the equations for spherical contact are dominated by the smaller radius object, this effect is somewhat less but still exists. Second, using the same constant coefficient for different loading histories is a approximation. For a higher fidelity, these effects need to be taken into account.

Chapter 7
Code Speedup and Parallelization

Planned and existing pebble bed reactors can have on the order of 100,000 pebbles. For some simulations, these pebbles need to be followed for long time periods, which can require computing billions of time-steps. Multiplying the time-steps required by the number of pebbles being computed over leads to the conclusion that large numbers of computations are required. These computations should be as fast as possible, and should be as parallel as possible, so as to allow relevant calculations to be done in a reasonable amount of time. This chapter discusses the process of speeding up the code and parallelizing it.

The PEBBLES program has three major portions of calculation. The first is determining which pebbles are in contact with other pebbles. The second computational part is determining the time derivatives for all the vectors for all the pebbles. The third computational part is using the derivatives to update the values. Overall, for calculation of a single time-step, the algorithm’s computation time is linearly proportional to the number of pebbles, that is O(n)¹ .

7.1 General Information about profiling

There are four different generic parts of the complete calculation that need to considered for determining the overall speed. The first consideration is the time to compute arithmetic operations. Modern processors can complete arithmetic operations in nanoseconds or fractions of nanoseconds. In the PEBBLES code, the amount of time spent on arithmetic is practically undetectable in wall clock changes. The second consideration is the time required for reading memory and writing memory. For main memory accesses, this takes hundreds of CPU clock cycles, so these times are on the order of fractions of microseconds (Drepper, 2007). Because of the time required to access main memory, all modern CPUs have on-chip caches, that contain a copy of the recently used data. If the memory access is in the CPU’s cache, the data can be retrieved and written in a small number of CPU cycles. Main memory writes are somewhat more expensive than main memory reads, since any copies of the memory that exist in other processor’s caches need to be updated or invalidated. So for a typical calculation like $a+b\to c$ the time spent doing the arithmetic is trivial compared to the time spent reading in $a$ and $b$ and writing out $c$.

The third consideration is the amount of time required for parallel programming constructs. Various parallel synchronization tools such as atomic operations, locks and critical sections take time. These take an amount of time on the same order of magnitude as memory writes. However, they typically need a read and then a write without any other processor being able to access that chunk of memory in between which requires additional overhead, and a possible wait if the memory address is being used by another process. Atomic operations on x86_64 architectures are faster than using locks, and locks are generally faster than using critical sections. The fourth consideration is network time. Sending and receiving a value can easily take over a millisecond for the round trip time. These four time consuming operations need to be considered when choosing algorithms and methods of calculation.

There are a variety of methods for profiling the computer code. The simplest method is to use the FORTRAN 95 intrinsics CPU_TIME and DATE_AND_TIME. The CPU_TIME subroutine returns a real number of seconds of CPU time. The DATE_AND_TIME subroutine returns the current wall clock time in the VALUES argument. With gfortran both these times are accurate to at least a millisecond. The difference between two different calls of these functions provide information on both the wall clock time and the CPU time between the calls. (For the DATE_AND_TIME subroutine, it is easiest if the days, hours, minutes, seconds and milliseconds are converted to a real seconds past some arbitrary time.) The time methods provide basic information and a good starting point for determining which parts of the program are consuming time. For more detailed profiling the oprofile (opr, 2009) program can be used on Linux. This program can provide data at the assembly language level which makes it possible to determine which part of a complex function is consuming the time. Non-assembly language profilers are difficult to accurately use on optimized code, and profiling non-optimized code is misrepresentative.

7.2 Overview of Parallel Architectures and Coding

Parallel computers can be arranged in a variety of ways. Because of the expense of linking shared memory to all processors, a common architecture is a cluster of nodes with each node having multiple processors. Each node is linked to other nodes via a fast network connection. The processors on a single node share memory. Figure 7.1 shows this arrangement. For this arrangement, the code can use both the OpenMP (Open Multi-Processing) (ope, 2008) and the MPI (Message Passing Interface) (mpi, 2009) libraries. MPI is a programming interface for transferring data across a network to other nodes. OpenMP is a shared memory programming interface. By using both programming interfaces high speed shared memory accesses can be used on memory shared on the node and the code can be parallelized across multiple nodes.

7.3 Lock-less Parallel O(N) Collision Detection

For any granular material simulation, which particles are in contact must be determined quickly and accurately for each time-step. This is called collision detection, though for pebble simulations it might be more accurately labeled contact detection. The simplest algorithm for collision detection is to iterate over all the other objects and compare each one to the current object for collision. To determine all the collisions using that method, $O\left(N^2\right)$ time is required.

An improved algorithm by Cohen et al. (1995) uses six sorted lists of the lower and upper bounds for each object. (There is one upper bound list and one lower bound list for each dimension.) With this algorithm, to determine the collisions for a given object, the bounds of the current objects are compared to bounds in the list—only objects that overlap the bounds in all three dimensions will potentially collide. This algorithm typically has approximately $O\left(Nlog\left(N\right)\right)$ time,² because of the sorting of the bounding lists (Cohen et al., 1995).

A third, faster method, grid collision detection, is available if the following requirements hold: (1) there is a maximum diameter of object, and no object exceeds this diameter, and (2) for a given volume, there is a reasonably small, finite, maximum number of objects that could ever be in that volume. These two constraints are easily satisfied by pebble bed simulations, since the pebbles are effectively the same size (small changes in diameter occur due to wear and thermal effects). A three-dimensional parallelepiped grid is used over the entire range in which the pebbles are simulated. The grid spacing $gs$ is set at the maximum diameter of any object (twice the maximum radius for spheres).

Two key variables are initialized, $grid\text{\_}count\left(x,y,z\right)$, the number of pebbles in grid locations x,y,z; and $grid\text{\_}ids\left(x,y,z,i\right)$, the pebble identification numbers ($ids$) for each x,y,z location. The $id$ is a unique number assigned to each pebble in the simulation. The spacing between successive grid indexes is $gs$, so the index of a given x location can be determined by $\lfloor \left(x-x_{min}\right)∕gs\rfloor$ where $x_{min}$ is the zero x index’s floor; similar formulas are used for y and z.

The grid is initialized by setting $grid\text{\_}count\left(:,:,:\right)=0$, and then the x,y,z indexes are determined for each pebble. The $grid\text{\_}count$ at that location is then atomically³ incremented by one and fetched. Because OpenMP 3.0 does not have a atomic add-and-fetch, the lock xadd x86_64 assembly language instruction is put in a function. The $grid\text{\_}count$ provides the fourth index into the $grid\text{\_}ids$ array, so the pebble $id$ can be stored into the $ids$ array. The amount of time to zero the $grid\text{\_}count$ array is a function of the volume of space, which is proportional to the number of pebbles. The initialization iteration over the pebbles can be done in parallel because of the use of an atomic add-and-fetch function. Updating the grid iterates over the entire list of pebbles so the full algorithm for updating the grid is $O\left(N\right)$ for the number of pebbles.

Once the grid is updated, the nearby pebbles can be quickly determined. Figure 7.2 illustrates the general process. First, index values are computed from the pebble and used to generate $xc$, $yc$, and $zc$. This finds the center grid location, which is shown as the bold box in the figure. Then, all the possible pebble collisions must have grid locations (that is, their centers are in the grid locations) in the dashed box, which can be found by iterating over the grid locations from $xc-1$ to $xc+1$ and repeating for the other two dimensions. There are $3^3$ grid locations to check, and the number of pebbles in them are bounded (maximum 8), so the time to do this is bounded. Since this search does not change any grid values, it can be done in parallel without any locks.

Therefore, because of the unique features of pebble bed pebbles simulation, a parallel lock-less $O\left(N\right)$ algorithm for determining the pebbles in contact can be created.

7.4 MPI Speedup

The PEBBLES code uses MPI to distribute the computational work across different nodes. The MPI/OpenMP hybrid parallelization splits the calculation of the derivatives and the new variables geometrically and passes the data at the geometry boundaries between nodes using messages. Each pebble has a primary node and may also have various boundary nodes. The pebble-primary-node is responsible for updating the pebble position, velocity, angular velocity, and slips. The pebble-primary-node also sends data about the pebble to any nodes that are the pebble boundary nodes and will transfer the pebble to a different node if the pebble crosses the geometric boundary of the node. Boundary pebbles are those close enough to a boundary that their data needs to be present in multiple nodes so that the node’s primary pebbles can be properly updated. Node 0 is the master node and does processing that is simplest to do on one node, such as writing restart data to disk and initializing the pebble data. The following steps are used for initializing the nodes and then transferring data between them:

1. Node 0 calculates or loads initial positions of pebbles.
2. Node 0 creates the initial domain to node mapping.
3. Node 0 sends domain to node mapping to other nodes.
4. Node 0 sends other nodes their needed pebble data.

Order of calculation and data transfers in main loop:

1. Calculate derivatives for node primary and boundary pebbles.
2. Apply derivatives to node primary pebble data.
3. For every primary pebble, check with the domain module to determine the current primary node and any boundary nodes.
   1. If the pebble now has a different primary node, add the pebble id to the transfer list to send to the new primary node.
   2. If the pebble has any boundary nodes, add the pebble id to the boundary send list to send it to the node for which it is a boundary.
4. If this is a time step where Node 0 needs all the pebble data (such as when restart data is being written), add all the primary pebbles to the Node 0 boundary send list.
5. Send the number of transfers and the number of boundary sends that this node has to all the other nodes using buffered sends.
6. Initialize three Boolean lists of other nodes that this node has:
   1. data-to-send-to with “true” if the number of transfers or boundary sends is nonzero, and “false” otherwise
   2. received-data-from to “false”
   3. received-the-number-of-transfers and the number-of-boundary-sends with “false.”
7. While this node has data to send to other nodes and other nodes have data to send to this node loop:
   1. Probe to see if any nodes that this node needs data from have data available.
      1. If yes, then receive the data and update pebble data and the Boolean lists as appropriate
   2. If there are any nodes that this node has data to send to, and this node has received the number of transfers and boundary sends from, then send the data to those nodes and update the Boolean data send list for those nodes.
8. Flush the network buffers so any remaining data gets received.
9. Node 0 calculates needed tallies.
10. If this is a time to rebalance the execution load:
    1. Send wall clock time spent computing since last rebalancing to node 0
    2. Node 0 uses information to adjust geometric boundaries to move work towards nodes with low computation time and away from nodes with high computation time
    3. Node 0 sends new boundary information to other nodes, and needed data to other nodes.
11. Continue to next time step and repeat this process until all time-steps have been run.

All the information and subroutines needed to calculate the primary and boundary nodes that a pebble belongs to are calculated and stored in a FORTRAN 95 module named network_domain_module. The module uses two derived types: network_domain_type and network_domain_location_type. Both types have no public components so the implementation of the domain calculation and the location information can be changed without changing anything but the module, and the internals of the module can be changed without changing the rest of the PEBBLES code. The location type stores the primary node and the boundary nodes of a pebble. The module contains subroutines for determining the location type of a pebble based on its position, primary and boundary nodes for a location type, and subroutines for initialization, load balancing, and transferring of domain information over the network.

The current method of dividing the nodes into geometric domains uses a list of boundaries between the z (axial) locations. This list is searched via binary search to find the nodes nearest to the pebble position, as well as those within the boundary layer distance above and below the zone interface in order to identify all the boundary nodes that participate in data transfers. The location type resulting from this is cached on a fine grid, and the cached value is returned when the location type data is needed. The module contains a subroutine that takes a work parameter (typically, the computation time of each of the nodes) and can redistribute the z boundaries up or down to shift work towards nodes that are taking less time computing their share of information. If needed in the future, the z-only method of dividing the geometry could be replaced with a full 3-D version by modifying the network domain module.

7.5 OpenMP Speedup

The PEBBLES code uses OpenMP to distribute the calculation over multiple processes on a single node. OpenMP allows directives to be given to the compiler that direct how portions of code are to be parallelized. This allows a single piece of code to be used for both the single processor version and the OpenMP version. The PEBBLES parallelization typically uses OpenMP directives to cause loops that iterate over all the pebbles to be run in parallel.

Some details need to be taken into consideration for the parallelization of the calculation of acceleration and torque. The physical accelerations imposed by the wall are treated in parallel, and there is no problem with writing over the data because each processor is assigned a portion of the total zone inventory of pebbles. For calculating the pebble-to-pebble forces, each processor is assigned a fraction of the pebbles, but there is a possibility of the force addition computation overwriting another calculation because the forces on a pair of pebbles are calculated and then the calculated force is added to the force on each pebble. In this case, it is possible for one processor to read the current force from memory and add the new force from the pebble pair while another processor is reading the current force from memory and adding its new force to that value; they could both then write back the values they have computed. This would be incorrect because each calculation has only added one of the new pebble pair forces. Instead, PEBBLES uses an OpenMP ATOMIC directive to force the addition to be performed atomically, thereby guaranteeing that the addition uses the latest value of the force sum and saves it before a different processor has a chance to read it.

For calculating the sum of the derivatives using Euler’s method, updating concurrently poses no problem because each individual pebble has derivatives calculated. The data structure for storing the pebble-to-pebble slips (sums of forces used to calculate static friction) is similar to the data structure used for the collision detection grid. A 2-D array exists where one index is the from-pebble and the other index is for storing $ids$ of the pebbles that have slip with the first pebble. A second array exists that contains the number of $ids$ stored, and that number is always added and fetched atomically, which allows the slip data to be updated by multiple processors at once. These combine to allow the program to run efficiently on shared memory architectures.

7.6 Checking the Parallelization

The parallelization of the algorithm is checked by running the test case with a short number of time steps (10 to 100). Various summary data are checked to make sure that they match the values computed with the single processor version and between different numbers of nodes and processors. For example, with the NGNP-600 model used in the results section, the average overlap of pebbles at the start of the run is 9.665281e-5 meters. The single processor average overlap at the end of the 100 time-step run is 9.693057e-5 meters, the 2 nodes average overlap is 9.693043e-5 meters, and the 12 node average overlap is 9.693029e-5 meters. The lower order numbers change from run to run. The start-of-run values match each other exactly, and the end-of-run values match the start of run values to two significant figures. However, the three different end-of-run values match to five significant digits. In short, the end values match each other more than they match the start values. The overlap is very sensitive to small changes in the calculation because it is a function of the difference between two positions. During coding, multiple defects were found and corrected by checking that the overlaps matched closely enough between the single processor calculation and the multiple processor calculations. The total energy or the linear energy or other computations can be used similarly since the lower significant digits also change frequently and are computed over all the pebbles.

7.7 Results

The data in Table 7.1 and Table 7.2 provide information on the time used with the current version of PEBBLES for running 80 simulation time steps on two models. The NGNP-600 model has 480,000 pebbles. The AVR model contains 100,000 pebbles. All times are reported in units of wall-clock seconds. The single processor NGNP-600 model took 251 seconds and the AVR single processor model took 48 seconds when running the current version. The timing runs were carried out on a cluster with two Intel Xeon X5355 2.66 GHz processors per node with a DDR 4X InfiniBand interconnect network. The nodes had 8 processors per node. The gfortran 4.3 compiler was used.

Significant speedups have resulted with both the OpenMP and MPI/OpenMP versions. A basic time step for the NGNP-600 model went from 3.138 seconds to 146 milliseconds when running on 64 processors. Since a full recirculation would take on the order of 1.6e9 time steps, the wall clock time for running a full recirculation simulation has gone from about 160 years to a little over 7 years. For smaller simulation tasks, such as simulating the motion of the pebbles in a pebble bed reactor during an earthquake, the times are more reasonable, taking about 5e5 time steps. Thus, for the NGNP-600 model, a full earthquake can be simulated in about 20 hours when using 64 processors. For the smaller AVR model, the basic time step takes about 34 milliseconds when using 64 processors. Since there are less pebbles to recirculate, a full recirculation would take on the order of 2.5e8 time steps, or about 98 days of wall clock time.

Chapter 8
Applications

The knowledge of the packing and flow patterns (and to a much lesser extent the position) of pebbles in the pebble bed reactor is an essential prerequisite for many in-core fuel cycle design activities as well as for safety assessment studies. Three applications have been done with the PEBBLES code. The major application is the computation of pebble positions during a simulated earthquake. Two other applications that have been done are calculation of space dependent Dancoff factors and calculation of the angle of repose for a HTR-10 simulation.

8.1 Applications in Support of Reactor Physics

8.1.1 Dancoff Factors

The calculation of Dancoff factors is an example application that needs accurate pebble position data. The Dancoff factor is used for adjusting the resonance escape probability for neutrons. There are two Dancoff factors that use pebble position data. The first is the inter-pebble Dancoff factor that is the probability that a neutron escaping from the fuel zone of a pebble crosses a fuel particle in another pebble. The second is the pebble-pebble Dancoff factor, which is the probability that a neutron escaping one fuel zone will enter another fuel zone without interacting with a moderator nuclide. Kloosterman and Ougouag (2005) use pebble location information to calculate the probability by ray tracing from fuel lumps until another is hit or the ray escapes the reactor. The PEBBLES code has been used for providing position information to J. L. Kloosterman and A. M. Ougouag’s PEBDAN program. This program calculate these factors as shown in Figure 8.2 which calculates them for the AVR reactor model.

8.1.2 Angle of Repose

The PEBBLES code was used for calculating the angle of repose for an analysis of the HTR-10 first criticality (Terry et al., 2006). The pebble bed code recirculated pebbles to determine the angle at which the pebbles would stack at the top of the reactor as shown in Figure 8.3, since this information was not provided, but was needed for the simulation of the reactor(Bäumer et al., 1990).

8.1.3 Pebble Ordering with Recirculation

During experimental work before the construction of the AVR, it was discovered that when the pebbles were recirculated, the ordering in the pebbles increased. Figures 8.4 and 8.5 show that this effect occurs in the PEBBLES simulation as well. The final AVR design incorporated indentations in the wall to prevent this from occurring.

8.2 Application to Earthquake modeling

The packing fraction of the pebbles in a pebble bed reactor can vary depending on the method of packing and the subsequent history of the packing. This packing fraction can affect the neutronics behavior of the reactor, since it translates into an effective fuel density. During normal operation, the packing fraction will vary only slowly, over the course of weeks and then stabilize. During an earthquake, the packing fraction can increase suddenly. This packing fraction change is a concern since packing fraction increase can increase the neutron multiplication and cause criticality concerns as shown by Ougouag and Terry (2001).

The PEBBLES code can simulate this increase and determine the rate of change and the expected final packing fraction, thus allowing the effect of an earthquake to be simulated.

8.2.1 Movement of Earthquakes

The movement of earthquakes has been well studied in the past. The magnitude of the motion of earthquakes is described by the Mercalli scale, which describes the maximum acceleration that a given earthquake will impart to structures. For a Mercalli X earthquake, the maximum acceleration is about 1 g. The more familiar Richter scale measures the total energy release of an earthquake (Lamarsh, 1983), which is not useful for determining the effect on a pebble bed core. For a given location, the soil properties can be measured, and using soil data and the motion that the bedrock will undergo, the motion on the surface can be simulated. The INL site had this information generated in order to determine the motion from the worst earthquake that could be expected over a 10,000 years period (Payne, 2003). This earthquake has roughly a Mercalli IX intensity. The data for such a 10,000 year earthquake are used for the simulation in this dissertation.

8.2.2 Method Of Simulation

The code simulates earthquakes by adding a displacement to the walls of the reactor. As well, the velocity of the walls needs to be calculated. The displacement in the simulation can be specified either as the sum of sine waves, or as a table of displacements that specifies the x, y, and z displacements for each time. At each time step both the displacement and the velocity of the displacement are calculated. When the displacement is calculated by a sum of sine functions, the current displacement is calculated by adding vector direction for each wave and the velocity is calculated from the sum of the first derivative of all the waves. When the displacement is calculated from a table of data, the current displacement is a linear interpolation of the two nearest data points in the table, and the velocity is the slope between them. The walls are then assigned the appropriate computed displacement and velocity. Figure 8.6 shows the total displacement for the INL earthquake simulation specifications that were used in this paper.

8.2.3 Earthquake Results

The results of two simulations carried out here show a substantially safer behavior than the Ougouag and Terry (2001) bounding calculations. The methodology was applied to a model of the PBMR-400 model and two different static friction coefficients were tested, 0.65 and 0.35. The packing fraction increased from 0.593 to 0.594 over the course of the earthquake with the 0.65 static friction model, with the fastest increase was from 0.59307 to 0.59356 and took place over 0.8 seconds. With the 0.35 static friction model, the overall increase was from 0.599 to 0.601. The fasted increase was from 0.59964 to 0.60011 in 0.8 seconds. This is remarkably small when compared to the bounding calculation packing fraction increase rate of 0.129 sec${}^{-1}$ in free fall.¹ Both computed increases and packing fraction change rates are substantially below the free fall bounding rate and packing fraction change of a transition from 0.60 to 0.64 in 0.31 seconds. The computed rate and the total packing fraction increase are in the range that can be offset by thermal feedback effects for uranium fueled reactors.

During the course of the earthquake, the boundary density fluctuations (that is the oscillations in packing fraction near a boundary) are observed to increase in amplitude. Figure 8.9 shows the packing fraction before the earthquake and after the earthquake in the radial direction. These were taken from 4 to 8 meters above the fuel outlet chute in the PBMR-400 model. All the radial locations have increased packing compared to the packing fraction before the earthquake, but the points that are at boundary density fluctuation peaks increase the most. This effect can be seen in figure 8.10, which shows the increase in packing fraction before the earthquake and after

A previous version of the positional data from the earthquake simulation was provided to J. Ortensi. This data was used by him to simulate the effects of an earthquake on a pebble bed reactor(Ortensi, 2009). Essentially, two factors cause an increase in reactivity. The first is the increased density of the pebbles and the second is due to the control rods being at the top of reactor, so when the top pebbles move down the control rod worth (effect) decreases. However, the reactivity increase causes the fuel temperature to rise, which causes Doppler broadening and more neutrons are absorbed by the ${}^{238}$U, which causes the reactivity to fall. Figure 8.11 shows an example of this.

8.2.4 Earthquake Equations

For each time-step, the simulation calculates both a displacement and a wall velocity.

For the sum of waves method, the displacement is calculated by:

where $t$ is the current time, $S$ is the time the wave starts, $p$ is the period of the wave, $c$ is the initial cycle of the wave, $o$ is the offset, and $D$ is the maximum displacement vector.

The velocity is calculated by:

For the tabular data, the displacement and velocity are calculated by:

where $T_i$ is the displacement at the $i$th time-step, $o$ is a number between 0 and 1 that specifies where 0 is the start of the time-step and 1 is the end, and $\delta$ is the time in seconds between time-steps.

With these displacements, the code then uses:

as the adjusted position and velocity.

Chapter 9
Construction of a Dust Production Framework

With the creation of the PEBBLES simulation, one issue that was examined was using the simulation to attempt to predict the volume of dust that would be produced by an operating pebble bed reactor. This is an important issue that could affect the choice of a pebble bed reactor versus a prismatic reactor for process heat applications. However, as this chapter and Appendix B will discuss, while the PEBBLES code has the force and motion data required for this simulation, the coefficients that would allow this information to be used have not been sufficiently robustly experimentally determined yet.

With the data provided by PEBBLES, equations to link the dust production to PEBBLES calculated quantities were examined. As shown in equation B.1, the volume of dust produced can be approximated if the normal force of contact, the length slide and the wear coefficients are known. The force of contact and the length slide are calculated as part of the PEBBLES simulation, so this method was used to calculate dust production for the AVR reactor. This resulted in an estimate of four grams of graphite dust produced per year as compared to the measured value of three kilograms of dust produced per year. Several possible causes of this were identified in the paper documenting this work (Cogliati and Ougouag, 2008). A key first issue as described by this dissertation’s previous work section is that there are no good measurements of graphite wear rates in pebble bed reactor relevant conditions (especially for a reactor run at AVR temperatures). A second issue is that the previous model of AVR was missing features including the reactor control rod nose cones and wall indentations. A third issue, identified after the paper’s publication, is that significant portions of the length traveled were due not to motion down through the reactor. Instead, much of the length that was tallied was due to pebbles vibrating. In the model used in the paper, this problem was about four times more severe than the current model, due to the new addition of slip correction via equation 5.2.

As an illustration of the general framework for dust production, a simple cylindrical vat the size of AVR was simulated. In this model, the outlet chute starts shrinking at height zero. The length (L), and length times normal force (NL) were tallied on intervals of 6 cm and the figures 9.1 to 9.5 show the results after recirculating 400 pebbles. The results are on a per pebble-pass basis. Two sets of static friction and kinetic friction pairs are used, one with a static friction coefficient of 0.35 and kinetic of 0.25, and the other with static of 0.65 and kinetic of 0.4. Figure 9.1 shows the calculated pebble-to-pebble lengths. Notice that for the 0.65 static friction simulation, about 10 meters of length traveled is occurring at the peak in a 6 cm long tally. Since the pebble is traveling about 0.06 m and has at most 12 pebble to pebble contacts, essentially all but a small portion of this length traveled is due to when the pebbles vibrate relative to each other. This vibration is caused by the impact of the pebbles coming from the inlet chutes and hitting the top of the bed. This likely is a true physical effect, which has not been discussed in literature this author is aware of. However, in order to obtain the correct magnitude of this vibrational effect, two things must be correct. First, the simulation must dissipate the vibration at the correct rate, and second, the pebbles must impact the bed at the correct velocity. Quality estimates of both will need to be made to finish the dust production work. Note that with the lower static friction value, the effect is even more pronounced.

Figure 9.2 is also expected to have a vibrational component, since only a small portion of pebbles should contact the wall, and therefore for a 6 cm tally, the length traveled by the average pebble should be much lower than 6 cm. As the chute is entered, the distance the average pebble travels increases in the 0.65 static friction case. Figure 9.3 shows the average normal contact force. The peak value for the pebble to surface values is due to the base of the ‘arch’ formed by pebbles. The curves for both the 0.65 static friction coefficient and the 0.35 static friction coefficient are approximately the same because the static friction force is not reaching the full Coulomb limit, so both have the same effective $\mu$.

Figure 9.4 shows the normal times force sums. For the 0.65 case, the peak is due to the vibrational impact. For the 0.35 case, the vibration travels deep into the reactor bed, producing dust throughout the reactor. Figure 9.5 shows the peak dust production coming in the base of the reactor, where the forces are the highest, and the greatest lengths are traveled next to the wall.

Chapter 10
Future Work

The dust production simulation requires both proper dust production wear coefficients, and properly determining the correct method of dealing with vibrational issues. It would be useful to determine the number of pebbles that need to be simulated to provide a correct representation of a full NGNP-600 sized reactor. Since the middle portions are geometrically similar, determining the amount of recirculation that is required to reach a geometrically asymptotic state might allow only a portion of the recirculation to be done. Those two changes might allow quicker simulation of full sized reactors. Finally, in order to allow sufficiently fast simulations on today’s computer hardware many approximations to the true behavior are done. In the future, some of these approximations maybe relaxed.

Chapter 11
Summary and Conclusions

Research results presented in this dissertation demonstrates a distinct element method that provides high fidelity and yet has reasonable run-times for many pebble fuel element flow simulations. The new static friction test will be useful for evaluating any implementation of static friction for spheres. The PEBBLES code produced for this dissertation has been able to provide data for multiple applications including Dancoff factor calculation, neutronics simulation and earthquake simulation. The new static friction model provides expected static friction behavior in the reactor including partial matching of the Janssen model predictions and correctly matching stability behavior in a pyramid. The groundwork has been created for predicting the dust production from wear in a pebble bed reactor once further experimental data is available. Future work includes potentially relaxing some of the physical approximations made for speed purposes when faster computing hardware exists, and investigating new methods for allowing faster simulations. This dissertation has provided significant enhancements in simulation of the mechanical movement of pebbles in a pebble bed reactor.

Bibliography

Appendix A
Calculation of Packing Fractions

For determining the volume of a sphere that is inside a vertical slice, the following formula can be used