Non-classical transport with

angular-dependent path-length distributions.

2: Application to pebble bed reactor cores

Richard Vasques

RWTH University, Aachen, Germany

Department of Mathematics

Center for Computational Engineering Science


Edward W. Larsen

University of Michigan, Ann Arbor, U.S.A.

Department of Nuclear Engineering and Radiological Sciences


We describe an analysis of neutron transport in the interior of model pebble bed reactor (PBR) cores, considering both crystal and random pebble arrangements. Monte Carlo codes were developed for (i) generating random realizations of the model PBR core, and (ii) performing neutron transport inside the crystal and random heterogeneous cores; numerical results are presented for two different choices of material parameters. These numerical results are used to investigate the anisotropic behavior of neutrons in each case and to assess the accuracy of estimates for the diffusion coefficients obtained with the diffusion approximations of different models: the atomic mix model, the Behrens correction, the Lieberoth correction, the generalized linear Boltzmann equation (GLBE), and the new GLBE with angular-dependent path-length distributions. This new theory utilizes a non-classical form of the Boltzmann equation in which the locations of the scattering centers in the system are correlated and the distance-to-collision is not exponentially distributed; this leads to an anisotropic diffusion equation. We show that the results predicted using the new GLBE theory are extremely accurate, correctly identifying the anisotropic diffusion in each case and greatly outperforming the other models for the case of random systems.


The pebble bed reactor (PBR), a concept which originated in Germany in the 1950’s, is a graphite-moderated, helium-cooled, very-high-temperature (generation IV) reactor. Several countries have addressed different possible PBR designs, amongst which we mention the HTR-PM [2] in China (following the successful test reactor HTR-10 [3]), the PBMR [4] in South Africa, and the MPBR [5] in the United States.

The fundamental PBR design is based on the use of spherical, same-sized fuel elements called pebbles. Each fuel pebble is made of pyrolytic graphite (the moderator), containing 10,000 microscopic fuel Tristructural-Isotropic (TRISO) particles. To achieve the desired reactivity, thousands of pebbles are piled on top of one another in a“random” manner (influenced by gravity) inside the cylindrical reactor core. They are dropped on top of this piling from charging tubes located at the top of the core, and move downward through the system in a dense granular flow. Due to this dynamic structuring, the exact locations of the pebbles inside the core at any given time are unknown.

Typically, the neutronic modeling of the geometrically random core is done by: (i) developing self-shielded multigroup cross sections for the pebbles, (ii) volume-averaging these cross sections over the entire core, including the helium-filled region between the pebbles (the atomic mix approximation), and (iii) introducing the spatially-homogenized cross sections into a diffusion code. This procedure leads to two concerns, as explicited next.

First, in the classic atomic mix model, the cross sections for a random heterogeneous medium are approximated by volume-averaging over the constituent materials, weighted by their respective volume fractions. This approximation is known [6, 7] to be accurate only when the chunk sizes of the constituent materials are small compared to a mean free path; however, the pebbles are O(1) mean free paths thick. In fact, it has been observed that neutron streaming in this type of system is strongly affected by the void spaces; to account for this effect, experimental and mathematical approaches were used to develop corrections for the diffusion coefficients obtained with atomic mix [8, 9]. Hence, the validity of the atomic mix approximation is called into question.

The second concern is related but subtly different: in a PBR core, does gravity affect the distance-to-collision in a direction-dependent (anisotropic) manner? In other words, the force of gravity (let us say it acts in the negative direction) causes the pebbles to arrange themselves in a certain manner, which is affected by the direction of this force. If one considers a typical arrangement of pebbles in a PBR core, the question is whether the chord length probability distribution function in the direction is different than in the -plane.

Currently, PBR cores are modeled using a diffusion approximation with an isotropic diffusion tensor, in which neutrons diffuse equally in all spatial directions. Previous research [10, 30] has indicated that, for different crystal arrangements of a PBR core, anisotropic diffusion effects can be found. In this work, we address this anisotropic behavior and also investigate its existence in PBR random structures.

The basic idea consists of using an angular-dependent, non-exponential ensemble-averaged probability distribution function to replace the true probability distribution function for the distance-to-collision. This leads to the new generalized linear Boltzmann equation (GLBE) derived in [11], and to its asymptotic diffusion limit: a diffusion equation with anisotropic diffusion coefficients. To investigate the accuracy of this extended theory, we have performed numerical simulations in model PBR cores, and compared the diffusion coefficients obtained numerically with those estimated by the atomic mix model (and its corrections) and by the new generalized theory. Overall, we find that the new theory yields extremely accurate results, correctly predicting anisotropic diffusion in each case and, more importantly, greatly outperforming the other models for the problems in random systems.

A summary of the remainder of the paper follows. In Section 2 we present the different formulations that model neutron transport in the PBR systems considered in this paper, and discuss how to obtain the diffusion coefficients. In Section 3 we show how the crystal and random realizations of the PBR model core were constructed. In Section 4 we discuss the Monte Carlo algorithm used to model neutron transport, and present the numerical results. In Section 5 we compare the theoretically estimated diffusion coefficients with the ones obtained numerically; and in Section 6 we conclude with a discussion.


The goal of this paper is to assess the accuracy of the new generalized theory in predicting the diffusion of neutrons in the interior of a PBR system; that is, away from boundary effects and strong packing fluctuations. In this section, we shortly describe the different formulations that model neutron transport and diffusion inside such systems, and present the expressions to estimate the diffusion coefficients.

Let the position vector be described by the cartesian coordinates , and let us write as the polar angle measured with respect to the (vertical) -axis and as the corresponding azimuthal angle. Introducing , we can write the vector


Now, consider a binary system composed of solid fuel (spherical) pebbles (material 1) immersed in a void background (material 2). In this case, for , we write the parameters:


Finally, defining

the path-length traveled by the neutron since (2.3)
its previous interaction (birth or scattering),

we make the following assumptions:

  • The physical system is both infinite and statistically homogeneous.

  • The system has azimuthal symmetry. (The probability distribution function for distance-to-collision depends only upon the polar angle.)

  • Neutron transport is monoenergetic.

  • Neutron transport is driven by a specified point source located at the origin and isotropically emitting neutrons per second.

  • The neutron flux as .

  • The ensemble-averaged total cross section , defined as

    is known.

  • Scattering is isotropic.

Assumption follows directly from the fact that, for any given PBR system (crystal or random), a realization containing pebbles with coordinates and a realization containing pebbles with coordinates , where , have the same probability of occurring. Therefore, when investigating the behavior of neutrons born in a fuel pebble in the interior of the system, we can assume without loss of generality that the mean-squared displacements in the and directions are the same.

2.1 The Atomic Mix Model

Given a single realization of the system, we can write the packing fraction of material 1 as


Then, defining the operator as the ensemble average over all possible realizations of the system, the atomic mix parameters are given by


The steady-state atomic mix equation [6] for this system is


where . Defining , the asymptotic diffusion approximation for Eq. (2.6) is


where D is the atomic mix diffusion coefficient. Notice that this model assumes classical transport, in which the probability distribution function of the distance traveled between collisions is an exponential. In this case, the mean and mean-squared free paths are given respectively by and .

2.2 Corrections to the Atomic Mix Diffusion Coefficient

In 1949, Behrens investigated the increase in the migration length of neutrons in a reactor caused by the presence of “holes” in the reactor [8]. In that work, “holes” are primarily understood to be the coolant spaces, due to the low density of the substances used as coolants. He also noticed that a small anisotropic effect occurred depending on the shape of the holes. For the case of pebble beds in which , he proposed that the isotropic diffusion coefficient D be given by


where is the total cross section of the solid material (pebbles), is the radius of a pebble, is the hole/material volume ratio of the system, D is given by (2.7b), and q is the quotient of the mean-squared path-length through the hole divided by the square of the mean path-length through the hole, estimated by


Later work [12, 13] emphasized the general validity of these equations for pebble bed problems. In 1980, Lieberoth & Stojadinović [9] revisited this theory. They developed a mock-up model of a pebble bed using steel balls and measured the coordinates of 3,024 sphere centers so that Monte Carlo games for neutron diffusion could be established. Using these results (as well as Monte Carlo calculations for crystal structures), they proposed an improved expression for q:


which also more closely relates to the theoretical value q = 2 obtained for randomly overlapping spheres [14]. Moreover, under the assumption that no correlation exists between the passage lengths in the holes and in the balls, they developed the following formula for the diffusion lengths:


This correction is used when more accurate estimates of neutron streaming in pebble bed type reactors are required [15, 16].

2.3 The New GLBE

Let us assume that the probability that a neutron will experience a collision while traveling an incremental distance in a direction depends on and ; that is, . Following [11], the new GLBE for this system is given by


The asymptotic diffusion approximation for Eq. (2.12) in the case of azimuthal symmetry is written as


where D is the diffusion coefficient in the direction given by the generalized theory, is the mean-squared free path of a neutron in the direction , and is the mean free path of a neutron. It is clear from these equations that the diffusion coefficients given by the new GLBE can differ in vertical and horizontal directions, accounting for anisotropic effects that may be present in the problem.

In the case of being independent of , such that , the diffusion equation reduces to the one in [17]:


where D represents a non-classical isotropic diffusion coefficient.

2.4 Exact Expressions for and Mean-Squared Displacements

Consider the diffusion equation below, taking place in an infinite system with a point source at the origin:


Bearing in mind that and as , we can manipulate this equation to obtain exact formulas for and for the mean-squared displacement of neutrons.

Operating on Eq. (2.15) by


we obtain the exact expression


Bearing in mind that represents the rate at which path-length is generated by neutrons in about , we see that


Introducing this result into Eq. (2.17) we obtain the exact expression


Notice that for the classic case (in which ) this expression reduces to the classic expression .

Next, multiplying Eq. (2.15) by and operating on it by Eq. (2.16), we obtain:

Integrating the first term on this equation by parts we get

We can rewrite Eq. (2.20a) as


where represents the mean-squared displacement of a neutron in the direction. This argument yields similar results for the and directions, giving us an exact expression for the diffusion coefficient in direction in terms of the mean-squared displacement:



In this section we discuss the construction of crystal and random systems that model a realization of a PBR core. The cylindrical geometry of the actual core is not addressed, since we are interested in analyzing neutron diffusion in the interior of the system (away from the boundaries).

3.1 Crystal Structures

Figure 1: Arrangement of pebbles in a horizontal layer with a given distance

There are several possible ways of piling identical hard spheres in a crystal-like structure. In particular, it has been shown [18, 19, 20] that the maximum packing fraction of identical spheres in a container is given by the face-centered cubic arrangement (FCC), with packing fraction . For this reason, we have based the crystal structures in this work in the FCC arrangement.

Let d be the diameter of a fuel pebble, and be the fixed distance between pebbles in the same layer, as shown in Figure 1. We place the first layer (A) of pebbles at the bottom of the system, and lock them in place. We then proceed to fill the system in a face-centered fashion; that is, positioning the second (B) and third (C) layers and sequentially repeating this structure. A finite example of this type of piling is shown in Figure 2. The height of the layer can be defined directly from the previous layers by


with .

For , this packing method yields the classic FCC structure, with coordination number (number of spheres contacted by a given sphere) 12. For the cases with , however, the

Figure 2: FCC crystal structure () in a box with side d
Figure 3: Crystal packings with (left); (center); (right)

“cubic” feature of these face-centered systems is lost; the structure resembles that of a face-centered orthorhombic. Moreover, since we do not allow pebbles to overlap each other, must not exceed a maximum value if the crystal-like structure of the packing is to be maintained:


The packing structures generated with have coordination number 6. The packing generated with is not a rotation of the one generated with ; it is rather a geometrically different structure, with coordination number 8. A diagram of each case is given in Figure 3. The graph in Figure 4 shows the packing fraction as a function of (in increments of 0.025).

Figure 4: Packing fractions in crystal structures with different values of

3.2 Random Structures

As summarized in [21], exhaustive experimental work has lead to believe that randomly packed spheres of the same diameter cannot have a packing fraction larger than ; in fact, recent analytical work yields the same results [22]. It has also been argued, however, that random packing itself is not well-defined [23, 24].

For the problem of PBR cores, the most accepted average packing fraction ranges from 0.60 to 0.62, values that were experimentally validated in [25]. However, the probability of occurrence of any single packing structure is not quantified, and as pointed out in [26], there exists neither experimental evidence nor theoretical proof to support the assertion that other packing arrangements within the core are impossible. In fact, it has been shown [27] that changes in the friction coefficients have a significant influence on the structuring of the pebbles; under certain loading circumstances, packing fractions of 0.59 are possible.

To generate the random packings presented in this work, we have used a variation of the ballistic deposition method [28, 29] introduced in [30]. In our ballistic algorithm, each pebble is released at a random point above a cubic box. It then follows a steepest descent trajectory until it reaches a position that is stable under gravity, in which case it has its coordinates stored. Given a (incomplete) realization of the system, we randomly drop and store the coordinates of 20 different tentative pebbles; we then choose the one with the lowest -coordinate to be added to the system, and discard the 19 remaining ones. Once a pebble is added to the system, its position is locked; that is, the pebble is frozen in place. Rearrangement of pebbles and/or cascading events cannot happen. No velocity or friction coefficients are taken into account; the only restriction is that a pebble can never, at any point of its trajectory, overlap the limits of the box or another pebble. Once the tentative pebble with the lowest -coordinate is added to the system, the process is repeated;

Figure 5: Example of a pebble bed random structure in a box with side d

pebbles continue to be added until the box is full. An example of a random piling obtained with this procedure is shown in Figure 5.

We have developed 100 different random packings in a box with side . Assuming an imaginary box with side inside the system, such that its walls are a distance away from the walls of the box, we can define the packing fraction as the ratio between the total volume of pebbles inside this imaginary box (including partial pebbles) and the volume of the imaginary box. For the 100 simulated random packings in this work, we have found the average packing fraction in the interior of the system to be , with a standard deviation of 0.0012.


Unfortunately, due to the statistical nature of the heterogeneous medium and its effect on the transport of neutrons, there is generally no analytical expression that can be used to obtain the mean and mean-squared free paths in PBR random systems. Nevertheless, there is a logical and straightforward set of steps that can be followed in order to numerically estimate this quantity.

If we consider a particle P that is born (or scatters) at a random point inside the pebble , the total distance that this particle will travel inside the pebbles before experiencing a collision can be sampled from the exponential distribution


where is the total cross section of the pebbles.

Let be the linepath starting at point along which P travels, and let be the length of inside the pebble . If , P will leak out of the pebble before experiencing a collision. It will then travel

Figure 6: Linepath of a particle between collisions

a distance along in the vaccum before entering a new pebble . If , the particle will leak out of without experiencing a collision and will travel some distance along in the vacuum before entering another pebble (). Eventually, if the particle does not leak out of the system, there will be a pebble in which , meaning that P will experience a collision within (Figure 6). The distance traveled by this particle between birth and collision is given by


In a given realization of this system, we use the following procedure:

  • Choose a pebble in the system.

  • Randomly choose a point inside this pebble.

  • Using Eqs. (4.1) and (4.2), calculate and .

  • Repeat this process for a large number of particles.

We tally the results obtained with this process in different directions , and then divide and by the number of linepaths generated to obtain and , the mean and mean-squared distance-to-collision in the fixed direction .

Problem d d d
1 1.0 0.99 0.01 0.99 1/4
2 2.0 0.995 0.005 0.9975 1/4
Table 1: Parametes for fuel pebbles with diameter d

For the two sets of parameters considered in this work, we assumed the background material in which the pebbles were piled to be vacuum; the parameters used for the material of the pebbles are given in Table 1. The neutron histories within the systems are determined by a Monte Carlo transport code as follows:

  • A neutron is born at a random point inside the fuel pebble.

  • A random direction of flight is defined.

  • The distance that travels inside the pebbles is sampled from Eq. (4.1).

  • undergoes a collision and all relevant information is stored.

  • If is absorbed, the history of ends and the algorithm goes back to step I.

  • If is scattered, the algorithm goes back to step II.

In this algorithm, we make use of the azimuthal symmetry of the system to improve our results involving the -plane. In short, this means that we obtain (without loss of generality) the result .

Problem 1 Problem 2
d d d d d d d d d
0.000 1.3501 3.7634 125.15 124.96 0.6751 0.9728 128.90 128.84
0.025 1.4003 4.0685 135.27 135.59 0.7003 1.0570 139.55 140.77
0.050 1.4497 4.3830 145.44 146.05 0.7250 1.1452 150.68 153.68
0.075 1.4982 4.7070 156.12 157.30 0.7492 1.2366 162.54 166.08
0.100 1.5452 5.0358 167.00 168.61 0.7728 1.3308 174.15 178.79
0.125 1.5909 5.3687 177.51 180.10 0.7956 1.4268 186.53 192.58
0.150 1.6349 5.7031 188.50 192.02 0.8176 1.5240 198.63 205.33
0.175 1.6769 6.0349 199.07 202.45 0.8387 1.6214 211.23 218.42
0.200 1.7165 6.3604 209.74 214.16 0.8585 1.7179 223.91 231.18
0.225 1.7538 6.6769 220.35 224.06 0.8770 1.8118 235.74 243.91
0.250 1.7879 6.9768 230.37 233.63 0.8942 1.9026 247.41 255.24
0.275 1.8189 7.2588 239.58 243.23 0.9096 1.9876 258.88 266.86
0.300 1.8459 7.5125 248.20 251.55 0.9233 2.0655 268.49 276.59
0.325 1.8690 7.7347 255.47 258.64 0.9348 2.1341 278.33 284.99
0.350 1.8873 7.9166 261.44 263.55 0.9439 2.1904 286.02 290.14
0.375 1.9004 8.0490 266.50 267.98 0.9504 2.2315 292.04 295.23
0.400 1.9074 8.1193 268.96 269.23 0.9540 2.2547 296.11 297.01
0.425 1.9080 8.1260 269.74 269.10 0.9542 2.2554 296.91 295.76
0.450 1.9007 8.0497 268.14 265.50 0.9507 2.2330 294.33 290.01
0.475 1.8851 7.8904 262.72 259.34 0.9428 2.1825 288.24 281.70
0.500 1.8596 7.6390 254.66 250.14 0.9300 2.1034 278.96 271.24
0.525 1.8225 7.2841 243.20 239.29 0.9115 1.9942 265.10 256.42
0.550 1.7723 6.8263 227.98 223.77 0.8864 1.8552 247.07 238.01
0.575 1.7063 6.2597 209.48 205.57 0.8533 1.6854 224.80 216.64
0.600 1.6209 5.5806 186.37 183.61 0.8106 1.4855 197.92 192.47
0.625 1.5110 4.7854 159.37 158.03 0.7557 1.2567 166.72 165.40
Table 2: Monte Carlo results in crystal structures with different values of

4.1 Numerical Results in Crystal Structures

For all crystal structures considered in this work, the model core consists of a periodically infinite structure of pebbles of diameter d. The histories of 3 neutrons were simulated in each system; Monte Carlo numerical results are given in Table 2. The statistical error (with 95% confidence) is less than 0.037 for all values of , , , and .

Figure 7: Ratio in crystal structures as a function of

As expected, we detect a clear anisotropic effect in each of these systems, as can be seen by the ratio depicted in Figure 7. This indicates that the diffusion coefficients must be different for the vertical and horizontal directions.

4.2 Numerical Results in Random Structures

For all random systems in this work, the packing of pebbles of diameter d took place in a cubic box with side d. We choose the pebble closest to the center of the packing structure to be the one where neutrons are born. Being interested in the behavior of neutrons in the interior of the system, we want to minimize the effect of the boundaries of the box (walls, top, bottom). According to the work in

Figure 8: Imaginary box positioned inside a random realization
Figure 9: Reinsertion of a particle in the imaginary box

[31], we need to consider pebbles that are three to five diameters off-walls in order to have a packing structure that is not influenced by the walls and by the bottom. For each random realization, we found that the fluctuations in packing fractions ceased being significant around a distance of two diameters from the boundaries, as was the case in the experiment performed in [9]. Nevertheless, according to the work in [31], one needs to consider pebbles that are three to five diameters off-walls in order to have a packing structure that is not influenced by the walls and by the bottom. Thus, in the Monte Carlo simulations performed, we allow neutrons to travel only inside the imaginary box with side d depicted in Figure 8.

The difference in the algorithm is in dealing with neutrons that leak out of this imaginary box. Let us assume that the center of the box is at the origin, and let us consider a particle that had its last collision at point inside a pebble A and that leaks out of the imaginary box through the plane d. First, defining the coordinates of the center of the pebble A as , we locate the pebble B with center at closest to the point . Then, we reinsert into the system at the point , as shown in Figure 9. Finally, we shift the whole system so that now the coordinates of the center of the box are , and proceed with the history of the particle. A similar process is used if the particle leaks through any of the other walls; we repeat this reinsertion and shifting procedure as many times as necessary. In other words, particles are traveling in an infinite “quasi-periodic” structure; here,

Problem 1 Problem 2
d d d d d d d d
Average 1.7053 6.2898 209.56 210.01 0.8527 1.7016 224.18 224.78
Error 0.095% 0.211% 0.201% 0.200% 0.080% 0.195% 0.172% 0.179%
Table 3: Monte Carlo results in random structures
Figure 10: Ratios in 100 random realizations

we use the term “quasi-periodic” because the system is not always shifted by the same values in a given direction.

We have simulated the histories of neutrons in each realization of the random system. The statistical error in each realization was found to be (with 97.5 confidence) less than 0.056 in Problem 1 and 0.032 in Problem 2 for all values of , , , and . The ensemble-averaged Monte Carlo results and the statistical errors (with 95 confidence) are given in Table 3. Each single realization presents a small anisotropic effect (see Figure 10). This anisotropy is not consistent: (i) we see that in about 25% of the realizations; and (ii) the ensemble-averaged values of and are less than 0.27% apart in both problems. These facts indicate that, on average, neutrons tend to travel further in the vertical direction than in the horizontal direction; however, this extra distance is very small.

Furthermore, we remark that we are working with an average packing fraction of 59.34%, which is about 2% smaller than the generally assumed average packing fraction in a PBR core (60%-62%). Thus, one expects the small anisotropic effects found in each realization to be even smaller in the interior of real PBR cores, since the pebbles are packed more closely and the void fraction is reduced.

It is not unreasonable to assume that, if enough realizations of the system are simulated, one should find that approaches a value of 1. This indicates that, at least for fuel pebbles in the interior of the system, it should not be necessary to worry about anisotropic diffusion. However, the diffusion of neutrons that are born in pebbles positioned close to the walls is anisotropic [32]; we do not address this issue here.

Monte Atomic Behrens Lieberoth Isotropic New
Carlo Mix Correction Correction GLBE GLBE
D Eq. (2.22) Eq. (2.7b) Eq. (2.8) Eq. (2.11) Eq. (2.14b) Eq. (2.13b)
D Eq. (2.22) Eq. (2.7b) Eq. (2.8) Eq. (2.11) Eq. (2.14b) Eq. (2.13c)
Table 4: Expressions to compute the diffusion coefficients for each model
Monte Atomic Mix “Old” New
Carlo and Corrections GLBE GLBE
0.000 0.4635 0.4628 0.4502 0.4705 0.4626 0.4646 0.4646 0.4646
0.025 0.4830 0.4841 0.4676 0.4903 0.4833 0.4842 0.4836 0.4855
0.050 0.5016 0.5037 0.4841 0.5093 0.5031 0.5039 0.5027 0.5063
0.075 0.5210 0.5250 0.5004 0.5282 0.5229 0.5236 0.5219 0.5271
0.100 0.5404 0.5456 0.5162 0.5467 0.5423 0.5432 0.5410 0.5476
0.125 0.5579 0.5660 0.5313 0.5646 0.5610 0.5624 0.5598 0.5677
0.150 0.5765 0.5873 0.5463 0.5825 0.5797 0.5814 0.5784 0.5875
0.175 0.5936 0.6036 0.5603 0.5993 0.5973 0.5998 0.5966 0.6063
0.200 0.6109 0.6238 0.5743 0.6160 0.6149 0.6176 0.6142 0.6244
0.225 0.6282 0.6388 0.5861 0.6304 0.6299 0.6345 0.6310 0.6416
0.250 0.6443 0.6534 0.5977 0.6444 0.6446 0.6504 0.6469 0.6573
0.275 0.6586 0.6686 0.6077 0.6567 0.6574 0.6651 0.6618 0.6719
0.300 0.6723 0.6814 0.6170 0.6680 0.6693 0.6783 0.6752 0.6845
0.325 0.6834 0.6919 0.6253 0.6782 0.6800 0.6897 0.6871 0.6950
0.350 0.6926 0.6982 0.6312 0.6854 0.6875 0.6991 0.6970 0.7034
0.375 0.7012 0.7051 0.6355 0.6907 0.6931 0.7059 0.7046 0.7086
0.400 0.7050 0.7058 0.6378 0.6935 0.6960 0.7095 0.7089 0.7105
0.425 0.7069 0.7052 0.6374 0.6930 0.6955 0.7098 0.7102 0.7091
0.450 0.7054 0.6984 0.6346 0.6896 0.6919 0.7059 0.7072 0.7031
0.475 0.6968 0.6879 0.6305 0.6845 0.6867 0.6976 0.6999 0.6931
0.500 0.6847 0.6725 0.6212 0.6732 0.6747 0.6846 0.6877 0.6785
0.525 0.6672 0.6565 0.6089 0.6581 0.6590 0.6661 0.6697 0.6589
0.550 0.6432 0.6313 0.5931 0.6389 0.6388 0.6419 0.6457 0.6344
0.575 0.6139 0.6024 0.5711 0.6122 0.6108 0.6114 0.6149 0.6045
0.600 0.5749 0.5664 0.5423 0.5776 0.5747 0.5738 0.5763 0.5688
0.625 0.5274 0.5229 0.5048 0.5334 0.5283 0.5278 0.5285 0.5266
Table 5: Diffusion coefficients in crystal structures: Problem 1.


In this section, we compare the diffusion coefficients obtained numerically with those estimated by the different models presented in Section 2.

Monte Atomic Mix “Old” New
Carlo and Corrections GLBE GLBE
0.000 0.2386 0.2385 0.2251 0.2455 0.2385 0.2402 0.2402 0.2402
0.025 0.2491 0.2513 0.2338 0.2565 0.2508 0.2516 0.2510 0.2528
0.050 0.2598 0.2650 0.2420 0.2672 0.2627 0.2633 0.2621 0.2656
0.075 0.2712 0.2771 0.2502 0.2780 0.2747 0.2751 0.2734 0.2785
0.100 0.2817 0.2892 0.2581 0.2886 0.2864 0.2870 0.2848 0.2914
0.125 0.2931 0.3025 0.2656 0.2989 0.2979 0.2989 0.2963 0.3041
0.150 0.3037 0.3139 0.2732 0.3093 0.3094 0.3107 0.3077 0.3165
0.175 0.3148 0.3256 0.2802 0.3191 0.3203 0.3222 0.3191 0.3284
0.200 0.3260 0.3366 0.2871 0.3289 0.3312 0.3335 0.3302 0.3401
0.225 0.3360 0.3476 0.2931 0.3373 0.3406 0.3443 0.3410 0.3509
0.250 0.3459 0.3568 0.2989 0.3456 0.3498 0.3546 0.3514 0.3611
0.275 0.3558 0.3667 0.3039 0.3528 0.3578 0.3642 0.3611 0.3703
0.300 0.3635 0.3745 0.3085 0.3595 0.3653 0.3729 0.3700 0.3785
0.325 0.3722 0.3811 0.3127 0.3655 0.3720 0.3805 0.3781 0.3853
0.350 0.3788 0.3842 0.3156 0.3698 0.3768 0.3868 0.3849 0.3905
0.375 0.3841 0.3883 0.3177 0.3729 0.3803 0.3913 0.3901 0.3938
0.400 0.3880 0.3892 0.3189 0.3746 0.3821 0.3939 0.3934 0.3949
0.425 0.3890 0.3875 0.3187 0.3743 0.3818 0.3940 0.3943 0.3933
0.450 0.3870 0.3813 0.3173 0.3723 0.3796 0.3915 0.3927 0.3889
0.475 0.3822 0.3735 0.3153 0.3693 0.3762 0.3858 0.3879 0.3816
0.500 0.3749 0.3646 0.3106 0.3625 0.3687 0.3769 0.3798 0.3713
0.525 0.3635 0.3516 0.3045 0.3537 0.3588 0.3646 0.3680 0.3579
0.550 0.3484 0.3356 0.2966 0.3423 0.3462 0.3488 0.3523 0.3418
0.575 0.3293 0.3174 0.2855 0.3266 0.3287 0.3292 0.3324 0.3227
0.600 0.3052 0.2968 0.2711 0.3065 0.3063 0.3054 0.3078 0.3008
0.625 0.2758 0.2736 0.2524 0.2810 0.2779 0.2772 0.2778 0.2759
Table 6: Diffusion coefficients in crystal structures: Problem 2.

Table 4 summarizes the expressions to compute the different diffusion coefficients. The Monte Carlo diffusion coefficients D and D are computed using the numerically calculated mean free path and mean-squared displacements and . The values obtained with the generalized theory (D, D, and D) use the numerical results for , , and . The atomic mix estimates (D, D, and D) use the Monte Carlo ensembe-averaged packing fraction for the random structures and the packing fractions depicted in Figure 4 for the crystal structures.

To investigate the accuracy of the theoretical estimates, we define the relative differences (errors) between the theoretical and Monte Carlo diffusion coefficients in the direction as


where can be replaced by and ; mc stands for the results obtained by monte carlo; and model represents the model being compared.

5.1 Estimates for Crystal Structures

The estimates for the diffusion coefficients in crystal structures are given in Tables 5 (Problem 1) and 6 (Problem 2), for different values of considered.

Figure 11: Percent relative errors in crystal structures for different values of

The results obtained with classic atomic mix consistently underestimate the diffusion coefficients in both directions, differing from the numerical results by large amounts: up to 10% in Problem 1 and 18% in Problem 2. On the other hand, the estimates computed with the other models yield much smaller differences.

The relative errors obtained by each model [as defined in Eq. (5.1)] are plotted in Figure 11. The Behrens correction presents the largest errors, despite being relatively accurate. The correction suggested by Lieberoth and the results obtained with the Isotropic GLBE perform similarly. Since DD for almost every choice of , the accuracy of both methods alternates accordingly to the anisotropy encountered in each case. Once again, we point out that these methods give isotropic diffusion coefficients, which do not model the anisotropic behavior of neutrons found in these systems.

However, the new GLBE correctly identifies the anisotropic behavior in every case; that is, it yields DD when DD, and DD when DD. It is an improvement over the other methods, with a maximum error of 0.89% in Problem 1 (against 2.24% for D, 1.90% for D, and 1.80% for D), and 2.18% in Problem 2 (against 4.09% for D, 3.57% for D, and 3.93% for D).

5.2 Estimates for Random Structures

The estimates for the diffusion coefficients in random structures are given in Table 7, as well as the relative errors of the theoretical estimates compared to

Monte Atomic Mix “Old” New
Carlo and Corrections GLBE GLBE
Problem D D D D D D D D
Coeff. 0.6144 0.6157 0.5617 0.6009 0.5990 0.6147 0.6146 0.6150
1 error
(%) - - 8.580 2.201 2.506 0.049 0.029 -
(%) - - 8.776 2.411 2.716 0.166 - 0.126
Coeff. 0.3286 0.3295 0.2809 0.3200 0.3214 0.3326 0.3324 0.3329
2 error
(%) - - 14.542 2.617 2.214 1.204 1.154 -
(%) - - 14.771 2.877 2.475 0.934 - 1.034
Table 7: Diffusion coefficients in random structures

the ensemble-averaged Monte Carlo results. Similarly to the results in crystal structures, classic atomic mix performs poorly. The estimates computed by the Behrens and Lieberoth corrections underestimate the diffusion coefficients, with relative errors ranging from 2.2% to 2.9%. More importantly, the estimates computed with the GLBE are a clear improvement over the other methods, presenting an excellent level of accuracy. They outperform atomic mix and its corrections by one order of magnitude in Problem 1 and by a factor of 2 in Problem 2.

As with the crystal structures, the new GLBE correctly identifies the anisotropic behavior in both random problems, yielding a larger diffusion coefficient in the vertical direction. It is slightly more accurate than the Isotropic GLBE in Problem 1, and as accurate as the Isotropic GLBE in Problem 2.


In this work, we have investigated anisotropic diffusion of neutrons in the interior of model pebble bed reactor (PBR) cores, in which pebbles are arranged in crystal or random structures. Using Monte Carlo simulations, we have found clear anisotropic behavior in the crystal structures; neutrons travel longer in certain directions, according to the geometry and the packing fraction of the system. We have also found a very small anisotropic behavior in the ensemble-averaged random structures.

We have used the diffusion approximation of a new generalized linear Boltzmann equation (new GLBE) to estimate anisotropic diffusion coefficients that can capture this anisotropic behavior. This new theory utilizes a non-classical form of the Boltzmann equation in which the locations of the scattered centers (pebbles) are correlated and the distance-to-collision is not given by an exponential. In order to successfully apply this theory, we need to estimate the mean and the (angular-dependent) mean-squared free paths of neutrons; we do this numerically. This extra information is all microscopic in nature; it is not a closure relation. We have shown that, at least for problems of the pebble bed kind, this new approach can correctly identify even very small anisotropic diffusion, which is a feature not encountered in any of the standard approaches currently in use.

We compare the accuracy of the new GLBE estimates for the diffusion coefficients against other methods: the atomic mix model and two of its diffusion corrections, and the standard GLBE. These methods yield isotropic diffusion coefficients; that is, they cannot capture any anisotropic behavior present in the systems. Nevertheless, the results obtained with the standard GLBE and with both corrections to atomic mix present a good level of accuracy.

We find that the new GLBE is generally an improvement over the other methods in estimating the diffusion coefficients in crystal structures. It correctly predicts the anisotropic behavior of every system, and its maximum relative errors are about half as big as the ones obtained with the other methods.

More importantly, for the random structures, the new and standard GLBE models were shown to be a great improvement over the other methods, with relative errors one order of magnitude smaller in one problem and about 50% smaller in the other. The new GLBE slightly outperforms the standard GLBE in one problem, and has similar accuracy in the other. However, once again, the new GLBE is capable of detecting even the very small anisotropic diffusion found in these random systems, which the standard GLBE cannot do.

The GLBE is more costly to simulate than the atomic mix approximation. Atomic mix only requires the cross sections of the constituent materials and their volume fractions to be known. The GLBE requires much more detailed information, which must be obtained by constructing realizations of the random system and developing an accurate estimate of the ensemble-averaged distribution function for distance-to-collision. However, the GLBE preserves important statistical properties of the original random system, such as the ensemble-averaged probability distribution function of the distance-to-collision. The new GLBE preserves the more general, angular-dependent version of these quantities, which makes it a systematically more accurate alternative to the current methods.

We have also shown that, in the case of diffusion in the interior of random PBR systems, one should not have to worry about anisotropic diffusion. The anisotropic effect was found to be very small for the tested packing fractions, and it is likely to be even smaller for higher packing fractions encountered in actual cores. Nevertheless, diffusion is anisotropic for neutrons that are born in pebbles close to the boundaries of the core (neutrons tend to travel longer distances in directions parallel to the boundary), and the new GLBE theory can be used to capture this behavior.

In future work, we intend to extend these results for problems with anisotropic scattering and energy-dependence. Although it is computationally expensive, an algorithm to obtain is straightforward; we intend to use it to generate results that can be applied to the new GLBE in order to estimate the criticality eigenvalue . Monte Carlo benchmark results for in the heterogeneous core will need to be developed; we also intend to do this.


  • [1]
  • [2] Z. Zhang, Z. Wu, Y. Xu, Y. Sun, and F. Li “Design of Chinese Modular High-Temperature Gas-Cooled Reactor HTR-PM.”2 International Topical Meeting on High Temperature Reactor Technology. Beijing, China; 2004.
  • [3] Z. Wu, D. Lin, and D. Zhong “The Design Features of the HTR-10.” Nucl. Eng. Design 218, 25 (2002).
  • [4] A. Koster, H.D. Matzner, and D.R. Nicholsi “PBMR Design for the Future.” Nucl. Eng. Design 222, 231 (2003).
  • [5] A.C. Kadak “MIT Pebble-Bed Reactor Project.” Nucl. Eng. Tech. 39, 95 (2007).
  • [6] E.W. Larsen, R. Vasques, and M.T. Vilhena “Particle Transport in the 1-D Diffusive Atomic Mix Limit.” International Topical Meeting on Mathematics and Computation, Supercomputing, Reactor Physics and Nuclear and Biological Applications. Avignon, France; 2005.
  • [7] R. Vasques “A Review of Particle Transport Theory in a Binary Stochastic Medium.” M.Sc. thesis. PPGMAp - Universidade Federal do Rio Grande do Sul (2005).
  • [8] D.J. Behrens “The Effect of Holes in a Reacting Material on the Passage of Neutrons.” Proc. Phys. Soc. A 62, 607 (1949).
  • [9] J. Lieberoth and A. Stojadinović “Neutron Streaming in Pebble Beds.” Nucl. Sci. Eng. 76, 336 (1980).
  • [10] D. Mathews, V. Thibulsky, and R. Chawla “Anisotropic Diffusion Effects in Deterministic Pebble-Bed Lattices.” Trans. Am. Nucl. Soc. 68, 20 (1993).
  • [11] R. Vasques and E.W. Larsen “Non-classical transport with angular-dependent path length distributions. 1: Theory.” arXiv:1309.4817 [math-ph] (2013).
  • [12] R.D. Neef “Ausbreitung von Neutronenwellen in statistischen Schüttungen von kugelförmigen Elementen.” JÜL-1273, Kernforschungszentrum, Jülich (1974).
  • [13] W. Scherer, H. Gerwin, and R.D. Neef “Theoretische Analyse des Kritischen HTR-Experiments KAHTER.” JÜL-1136-RG, Kernforschungszentrum, Jülich (1974).
  • [14] W.C. Strieder and S. Prager “Knudsen Flow through a Porous Medium.” Phys. Fluids 11, 2544 (1968).
  • [15] W. Bernnat and W. Feltes “Models for Reactor Physics Calculations for HTR Pebble-Bed modular Reactors.” Nucl. Eng. Design 222, 331 (2003).
  • [16] T. Williams, M. Rosselet, and W. Scherer (editors) “Critical Experiments and Reactor Physics Calculations for Low-Enriched High Temperature Gas Cooled Reactors.” IAEA-TECDOC-1249, Vienna, Austria (2001).
  • [17] E.W. Larsen and R. Vasques “A Generalized Linear Boltzmann Equation for Non-Classical Particle Transport.” J. Quant. Spectrosc. Radiat. Transfer 112, 619 (2011).
  • [18] T.C. Hales “A proof of the Kepler conjecture.” Ann. Math. (Second Series) 162, 1065 (2005).
  • [19] T.C. Hales and S.P. Ferguson “A formulation of the Kepler conjecture.” Discrete Comput. Geom. 36, 21 (2006).
  • [20] T.C. Hales and S.P. Ferguson The Kepler Conjecture: The Hales-Ferguson Proof. Springer, New York, 2011.
  • [21] From the molecular physics correspondent “What is Random packing.” Nature 239, 488 (1972).
  • [22] C. Song, P. Wang, and H.A. Makse “A Phase Diagram for Jammed Matter.” Nature 453, 629 (2008).
  • [23] S. Torquato, T.M. Truskett, and P.G. Debenedetti “Is Random Close Packing of Spheres Well Defined?” Phys. Rev. Lett. 84, 2064 (2000).
  • [24] S. Torquato Random Heterogeneous Materials: Microstructure and Macroscopic Properties. Springer-Verlag, New York, 2002.
  • [25] M-M. El-Wakil Nuclear Energy Conversion. International Textbook Company, American Nuclear Society, 1982.
  • [26] A.M. Ougouag and W.K. Terry “A Preliminary Study Of The Effect Of Shifts In Packing Fraction On K-Effective In Pebble-Bed Reactors.” International Meeting on Mathematical Methods for Nuclear Applications. Salt Lake City, United States of America; 2001.
  • [27] J.J. Cogliati and A.M. Ougouag “Pebbles: A Computer Code for Modeling Packing, Flow, and Re-Circulation of Pebbles in a Pebble Bed Reactor.” 3 International Topical Meeting on High Temperature Reactor Technology. Johannesburg, South Africa; 2006.
  • [28] R. Jullien and P. Meakin “Simple Three-Dimensional Models for Ballistic Deposition with Restructuring.” Europhys. Lett. 4, 1385 (1987).
  • [29] P. Meakin and R. Jullien “Simulation of Small Particle Penetration in a Random Medium.” J. Phys. 51, 2673 (1990).
  • [30] R. Vasques “Anisotropic Diffusion of Neutral Particles in Stochastic Media” Ph.D. thesis. Department of Mathematics - University of Michigan (2009).
  • [31] D. Bedenig “Probleme des Kugelhaufens und Körniger Schüttungen.” Proc. THTR Symposium. Jülich, Germany; 1968.
  • [32] R. Vasques “Estimating Anisotropic Diffusion of Neutrons Near the Boundary of a Pebble Bed Random System.” International Conference on Mathematics and Computational Methods Applied to Nuclear Science & Engineering. Sun Valley, United States of America; 2013.
  • [33]
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description