Nonclassical transport with
angulardependent pathlength distributions.
2: Application to pebble bed reactor cores
Richard Vasques
RWTH University, Aachen, Germany
Department of Mathematics
Center for Computational Engineering Science
vasques@mathcces.rwthaachen.de
and
Edward W. Larsen
University of Michigan, Ann Arbor, U.S.A.
Department of Nuclear Engineering and Radiological Sciences
edlarsen@umich.edu
ABSTRACT
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 angulardependent pathlength distributions. This new theory utilizes a nonclassical form of the Boltzmann equation in which the locations of the scattering centers in the system are correlated and the distancetocollision 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.
1. INTRODUCTION
The pebble bed reactor (PBR), a concept which originated in Germany in the 1950âs, is a graphitemoderated, heliumcooled, veryhightemperature (generation IV) reactor. Several countries have addressed different possible PBR designs, amongst which we mention the HTRPM [2] in China (following the successful test reactor HTR10 [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, samesized fuel elements called pebbles. Each fuel pebble is made of pyrolytic graphite (the moderator), containing 10,000 microscopic fuel TristructuralIsotropic (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 selfshielded multigroup cross sections for the pebbles, (ii) volumeaveraging these cross sections over the entire core, including the heliumfilled region between the pebbles (the atomic mix approximation), and (iii) introducing the spatiallyhomogenized 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 volumeaveraging 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 distancetocollision in a directiondependent (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 angulardependent, nonexponential ensembleaveraged probability distribution function to replace the true probability distribution function for the distancetocollision. 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.
2. FORMULATIONS
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
(2.1) 
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:
(2.2a)  
(2.2b) 
Finally, defining
the pathlength 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 distancetocollision 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 ensembleaveraged 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 meansquared 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
(2.4) 
Then, defining the operator as the ensemble average over all possible realizations of the system, the atomic mix parameters are given by
(2.5a)  
(2.5b) 
The steadystate atomic mix equation [6] for this system is
(2.6) 
where . Defining , the asymptotic diffusion approximation for Eq. (2.6) is
(2.7a)  
(2.7b) 
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 meansquared 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
(2.8) 
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 meansquared pathlength through the hole divided by the square of the mean pathlength through the hole, estimated by
(2.9) 
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 mockup 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:
(2.10) 
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:
(2.11)  
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
(2.12)  
The asymptotic diffusion approximation for Eq. (2.12) in the case of azimuthal symmetry is written as
(2.13a)  
(2.13b)  
(2.13c) 
where D is the diffusion coefficient in the direction given by the generalized theory, is the meansquared 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]:
(2.14a)  
(2.14b) 
where D represents a nonclassical isotropic diffusion coefficient.
2.4 Exact Expressions for and MeanSquared Displacements
Consider the diffusion equation below, taking place in an infinite system with a point source at the origin:
(2.15) 
Bearing in mind that and as , we can manipulate this equation to obtain exact formulas for and for the meansquared displacement of neutrons.
Operating on Eq. (2.15) by
(2.16) 
we obtain the exact expression
(2.17) 
Bearing in mind that represents the rate at which pathlength is generated by neutrons in about , we see that
(2.18)  
Introducing this result into Eq. (2.17) we obtain the exact expression
(2.19) 
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:
(2.20a)  
Integrating the first term on this equation by parts we get  
(2.20b) 
We can rewrite Eq. (2.20a) as
(2.21) 
where represents the meansquared 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 meansquared displacement:
(2.22) 
3. CONSTRUCTING PBR MODELS
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
There are several possible ways of piling identical hard spheres in a crystallike 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 facecentered 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 facecentered 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
(3.1) 
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
“cubic” feature of these facecentered systems is lost; the structure resembles that of a facecentered orthorhombic. Moreover, since we do not allow pebbles to overlap each other, must not exceed a maximum value if the crystallike structure of the packing is to be maintained:
(3.2) 
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).
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 welldefined [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;
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.
4. MONTE CARLO NUMERICAL RESULTS
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 meansquared 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
(4.1) 
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
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
(4.2) 
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.

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 meansquared distancetocollision 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 
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 
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 .
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
[31], we need to consider pebbles that are three to five diameters offwalls 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 offwalls 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 “quasiperiodic” structure; here,
Problem 1  Problem 2  

d  d  d  d  d  d  d  d  
Ensemble  
Average  1.7053  6.2898  209.56  210.01  0.8527  1.7016  224.18  224.78 
Statistical  
Error  0.095%  0.211%  0.201%  0.200%  0.080%  0.195%  0.172%  0.179% 
we use the term “quasiperiodic” 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 ensembleaveraged 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 ensembleaveraged 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) 
Monte  Atomic Mix  “Old”  New  

Carlo  and Corrections  GLBE  GLBE  
D  D  D  D  D  D  D  D  
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 
5. ESTIMATED DIFFUSION COEFFICIENTS
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  
D  D  D  D  D  D  D  D  
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 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 meansquared 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 ensembeaveraged 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
(5.1) 
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.
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  
Diffusion  
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    
error  
(%)      8.776  2.411  2.716  0.166    0.126  
Diffusion  
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    
error  
(%)      14.771  2.877  2.475  0.934    1.034 
the ensembleaveraged 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.
6. DISCUSSION
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 ensembleaveraged 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 nonclassical form of the Boltzmann equation in which the locations of the scattered centers (pebbles) are correlated and the distancetocollision is not given by an exponential. In order to successfully apply this theory, we need to estimate the mean and the (angulardependent) meansquared 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 ensembleaveraged distribution function for distancetocollision. However, the GLBE preserves important statistical properties of the original random system, such as the ensembleaveraged probability distribution function of the distancetocollision. The new GLBE preserves the more general, angulardependent 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 energydependence. 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.
References
 [1]
 [2] Z. Zhang, Z. Wu, Y. Xu, Y. Sun, and F. Li “Design of Chinese Modular HighTemperature GasCooled Reactor HTRPM.”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 HTR10.” 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 PebbleBed Reactor Project.” Nucl. Eng. Tech. 39, 95 (2007).
 [6] E.W. Larsen, R. Vasques, and M.T. Vilhena “Particle Transport in the 1D 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 PebbleBed Lattices.” Trans. Am. Nucl. Soc. 68, 20 (1993).
 [11] R. Vasques and E.W. Larsen “Nonclassical transport with angulardependent path length distributions. 1: Theory.” arXiv:1309.4817 [mathph] (2013).
 [12] R.D. Neef “Ausbreitung von Neutronenwellen in statistischen Schüttungen von kugelförmigen Elementen.” JÜL1273, Kernforschungszentrum, Jülich (1974).
 [13] W. Scherer, H. Gerwin, and R.D. Neef “Theoretische Analyse des Kritischen HTRExperiments KAHTER.” JÜL1136RG, 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 PebbleBed modular Reactors.” Nucl. Eng. Design 222, 331 (2003).
 [16] T. Williams, M. Rosselet, and W. Scherer (editors) “Critical Experiments and Reactor Physics Calculations for LowEnriched High Temperature Gas Cooled Reactors.” IAEATECDOC1249, Vienna, Austria (2001).
 [17] E.W. Larsen and R. Vasques “A Generalized Linear Boltzmann Equation for NonClassical 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 HalesFerguson 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. SpringerVerlag, New York, 2002.
 [25] MM. ElWakil 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 KEffective In PebbleBed 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 ReCirculation 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 ThreeDimensional 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]