AN IMPROVED DOUBLE HETEROGENEITY MODEL FOR PEBBLE-BED REACTORS IN DRAGON-5

In pebble-bed reactors, the fuel is contained in small grains, which are included in a graphite matrix. Some burnable poison particles may also be present. In this work, an additional ’double heterogeneity’ model is introduced in the DRAGON5 lattice code. The model is based on the legacy work of She (INET) and has been improved to overcome intrinsic limitations. It is based on a simplified physical model whereas the two already existing models were based on the collision probability analysis or on renewal theory. The advantage of this new model is its simplicity to implement. The theory shows that the correction suggested in the original model should not be arbitrary, but a constant equivalent particle fraction of 0.63. Numerical comparison between the models is generally good. However it does not support the theory of a constant equivalent fraction. Additional work is needed to reduce the discrepancy between the models in some cases.


CONTEXT AND INTRODUCTION
In pebble-bed reactors, the fuel is contained in small grains, which are included in a graphite matrix. Here, the fuel dispersion reactors is studied and a new 'double heterogeneity' model is introduced in the DRAGON5 lattice code, based on the legacy work of She (INET). With a deterministic approach, a direct and exact geometry representation would be too expensive, and a homogenized mixture would be highly inaccurate. Several models have been introduced in the past, one based on the collision probability analysis (Hébert [1]), one based on renewal theory (Sanchez and Pomraning [2]), and finally, one based on a simplified physical model (She [3]). The later is very easy to implement, but some of its extensions to resonance self-shielding and fuel depletion had to be tested. DRAGON5 code already included the first two models, the later was thus programmed as proposed by She. The advantage of using the same code is that no bias is introduced when comparisons are done. The only differences come from the calculations of the equivalent cross-sections in the composite material, not the resonance treatment or the tracking options for examples.
In this paper, the model of She is briefly presented in section 2. The computation of the source term needed to solve the transport equation and the flux is detailed both for single type and multiple types of particles in the composite. Then, in section 3 an improvement to a major limit of the model is proposed. Test cases results are presented in section 4, before conclusions are drawn.

THEORY
As mentioned previously, several approaches can be found in the literature to take into account properly the double heterogeneity of geometries such as in pebble bed. The main goal is to compute the average cross-sections in the composites and to be able to obtain the flux in each of their components. This section presents the theory and model proposed by She et al. [3] harmonized with the notation used by Hebert [4]. It also defines several quantities needed by the different steps of DRAGON calculations and for the improvement presented in Section 3. The model proposed by She is introduced in a subroutine of DRAGON which computes the equivalent cross-sections of a composite. This subroutine is called by the self-shielding modules as well as the flux solver.

Unique Type of Microstructure in a Composite
First, the case with only one type of microstructure in the composite is studied, and the basic assumptions are presented. Additionally to Hébert notation to remove all possibilities of mistakes, notations have the superscript '1' when referring to a unique microstructure composite. To compute the different properties of the composite, She et al. [3] make the assumption that the composite i behaves as a cylinder with a spherical grain in its centre ( Fig. 1 and 2 in Ref [3]). The neutrons travel through the cylinder parallel to its axe. Thus, they enter on one round side, travel through the matrix, then the grain of type j (all its layers) and the matrix again, to finally exit on the other round side. Some of then collide along the way, some not, which defines p 1 (i, j), the probability for a neutron to go through without collision. Then, the total equivalent macroscopic cross-sectioñ Σ 1 i in a composite i is given by: where the length of the cylinder, L 1 (i), is chosen to keep the volume ratio of microstructure / matrix. To simplify the notations, all the reference to the grain type j, are removed in this section since only one type of grain is considered. The conservation of reaction rates equation (Eq. 1 of [3]) is then given by: where the indexes i, dil and k represents respectively: the composite, its matrix and the layer of the microstructure j,φ is the average flux in the whole composite, φ dil (i) and φ k (i) the flux in the different subregions of the composite, V , V 1 dil and W k the total, matrix and microstructure subregion volumes respectively. This can be rewritten as follows: where f are the volume fractions of the different subvolumes, with f 1 (i) = k f 1 k (i), and S are the self-shielding factors of the matrix and the layers of the microstructure j. Moreover, the total equivalent macroscopic cross-sectionΣ i in a composite i can also be expressed as the average of the equivalent cross-section of all its subvolumes: A comparison of the Eq. 2 to 4 shows that the spatial self-shielding factors are: As She et al mentioned in [3], the self-shielding can be applied to all reaction types. Thus, using notations from [4], the within-group scattering equivalent macroscopic cross-sectionΣ w i in a composite i is given by: Using the definition of the spatial self-shielding factor in Eq. 5 to the within-group scattering equivalent microstructure cross-sectionsΣ w,k (i, j), and Eq. 24 and 25 of Ref. [3], we can write : where p 1 k (i, j) and p 1 dil (i) are referred as p k j in Ref. [3]. Similarly to p 1 i , p 1 k (i) is obtained by following the path the neutrons have to travel to have their first collision in layer k of microstructure j of composite i. Note that the meaning of indexes j and k are switched from Ref. [3]. Moreover, the conservation of neutrons in terms of collision probabilities can be rewritten as follows: Using Eq. 7 and 8 in Eq. 6, we obtain: By taking a closer look at the previous equation, we can see that the within-group scattering equivalent macroscopic cross-sections is equal to the total equivalent cross-sections time a factor. This factor is defined as a weighted average on all microstructures (and matrix) of the ratio withingroup scattering divided by total cross-section. The weight is the probability of first collision in each subregion. Note that the weights (p 1 k (i) and p 1 dil (i)) are normalized by their sum which is the total collision probability (1 − p 1 (i)).
By definition, the source is the product of a cross-section (scattering or fission) and a flux. When the equivalent cross-section for the composite is computed, the reaction rates are preserved, which means that the sources also are. This leads to the following simple equation for the equivalent source : Note that the source has to be computed with the local flux, NOT the average composite flux to keep the previous equation valid.
Once all the equivalent values are obtained for all composites, the average flux can be computed.
To recover the flux in each matrix and layers of the microstructures, Eq. 5 is used. Local fluxes can then be used to perform depletion calculations, for example. Note that if there is no microstructure (homogeneous macro-volume), then in Eq. 10 and 11, p 1 k (i) = 0 and f 1 k (i) = f (i) = 0. Thus, according to Eq. 9, p 1 dil (i) and 1 − p 1 (i) are equal, which leads to:

Several Types of Microstructures in a Composite
For a composite with several types of microstructures, She et al propose the following procedure. First, the main assumption is that the grains do not have shadowing effect on each other, regardless of their number and type. The matrix can then be split between them accordingly to their volume fraction f (i, j) in the composite i. This procedure defines a representative volume for each microstructure j with a total volume V eq (i, j), a microstructure fraction f eq (i, j), and its layers fractions f eq k (i, j) with: where f (i) is the volume fraction of all microstructures together, and W k (i, j) are the volumes of each layer in the microstructure. Since there is no shadowing effect according to the main assumption, the composite equivalent macroscopic cross-section is then given by the average of the equivalent macroscopic cross-section on all representative volumes V eq (i, j): Similarly, for the composite equivalent within-group macroscopic cross-section and sources (excluding within-group scattering), we have: It is highly important here to note that the equivalent valuesΣ 1 w,i (j),Σ 1 w,i (j) and q as * ,1 i (j) are computed with the equivalent fraction f (i) (Eq.12) and NOT their actual volume fraction. The procedure presented in the previous section is followed for each microstructure independently, before the average is made. Again, once all the equivalent values are obtained for all composites, the average flux can be computed. To recover the flux in each layers of the microstructures, Eq. 5 is used: For the matrix, a similar equation in the representative volume for each microstructure is used. This leads to a flux value potentially different for each microstructure; the average is then performed according to the volume fraction as follows: With 12, 13 Eq. 16 and 17 , the flux distribution can be recovered in all subregions (matrix and layers of microstructures) and used for further calculations such as reaction rates, depletion...

LIMITATIONS AND IMPROVEMENT
As mentioned by She et al [3], the equivalent volume approach using a cylinder has an intrinsic default. Indeed, no verification is done on the length vs. radius of the cylinder. As the authors underlined, this becomes an issue when the microstructure volume fraction becomes small, then the cylinder become too long, much longer than the composite itself. To overcome this issue, She et al propose to increase the radius of the microstructure (R ), as if a layer of matrix mixture was applied around it. The new volume fraction f of this augmented microstructure is then larger, which means that the length of its cylinder is reduced (L ). The relative dimensions between initial and augmented cylinders are presented in Ref. [3]. She et al showed on an example that this correction is mandatory when the fraction is under 5%. Similar results will be presented later. They also suggest a value between 5 and 10% for the minimum arbitrary volume fraction f . In the remaining of this section, we would like to explore more what a specific arbitrary choice actually means in terms of physics and approximations done. On the one hand, as mentioned before, when the volume fraction is small and the length of the cylinder is artificially very long, it means the neutron have to travel through a very long path through the matrix before and after they get to the microstructure. This leads to an excessive and not realistic absorption in the matrix. On the other hand, the length of the cylinder could be chosen as short as possible to fit the pseudo-microstructure in it, which is its diameter. In that case, the pseudo-volume fraction is 2/3. This approach is closer to Hebert or Sanchez where the authors approximate the representative volume as a sphere of matrix around spheric microstructures. Now, one could ask which approach does make more sense physically. The former may not physically be acceptable and the later has the shortest distance of matrix. To reduce the approximations done in the cylinder model, the better choice should reflect the actual length travelled by a neutron between two microstructures. Even though the general approximation is to suppose that all microstructures are independent and randomly dispersed in the matrix, it would be more realistic if the length of the cylinder is equal to the average distance between microstructure. Then, the relative radius of the microstructure and the equivalent cylinder would be better suited to represent the number of neutrons in the composite which travel across the microstructure or only through the matrix along the cylinder length.
The first step is to calculate the average distance between each microstructure of one type. For that, we assume momentarily that the microstructures are on a regular tetrahedron pattern, in which all the closest microstructures are at the same distance. In that case, when spheres touch each other, the maximum density is obtained: where L is the length of the tetrahedron edges and of the equivalent cylinder we are looking for, V S is the volume occupied by the sphere in the tetrahedron, V T is the tetrahedron volume. For the same edges length, but smaller spheres (radius R), the volume fraction is given by: Thus for the cylinder shape, the volume fraction f and the equivalent fraction f are given by: The previous equation shows that the 'best' equivalent volume fraction is actually a constant, close to the maximum.

Test Case Description
To test the implementation and the improvement of the double heterogeneity model presented in the previous sections, pebbles similar to those used in the HTR-10 reactor [6] are considered. They consist in a spherical composite of carbon matrix with particles. The composite is 5cm diameter covered by a layer of the same carbon matrix (0.5 cm thick). The particles include a core of U O 2 fuel with a diameter of 250µm. Similar cases with additional burnable poison particles of 10 B 4 C are also considered. Both particles content and geometry are described completely in Ref. [5].

Simulations
Computations were performed with DRAGON using a Jeff 3.2 library and SHEM361 energy group structure. The results for pebbles with a single type of particle are presented in Table 1. Several amounts of uranium particles (w in grams per ball) were simulated with the method of Hébert (HEBE), Sanchez (SAPO) and She (SLSI). For the SLSI method, several minimum equivalent fractions of particles were tested: none, 5%, 10%, 20%, 35% and 63%. The results show that Hébert and Sanchez methods give very similar results. The minimum fraction has a large influence on the results (several hundreds of pcm between 5% and 63%). Similar tests performed with XPZ with ENDF-B7.0 by the author of Ref. [5] show the same phenomenon. When compared to the results with two Monte-Carlo based codes (SERPENT with Jeff 3.1.1 and RMC with ENDF-B7.0), the theoretical fraction of 63% calculated previously do not always give the closest results. By interpolation of the results, an approximate equivalent fraction of 1.5 time the real fraction gives the same results as SERPENT. However, when compared to the XPZ results, an equivalent fraction between 5 and 10% is the closest.  Table 2 presents the results for the simulation with burnable absorbers: B 4 C particles * . In this case, RMC results presented in Ref [5] are taken as reference. The larger the amount of absorbent is, the larger the discrepancy becomes. The influence of the minimum equivalent fraction is even more important with these cases, where the equivalent theoretical fraction of 63% gives the closest results.

CONCLUSIONS
The method developed by She to simulate the double heterogeneity in composite, such as in pebblebed fuel, has been successfully implemented in DRAGON. Moreover, the limitations of the method regarding small fraction of particles in composites have been studied. In this paper, it has been demonstrated theoretically that the equivalent fraction used to define the representative volume is actually a constant close to the maximum that a cylinder approach can allow. The numerical results support the theory and the improvement proposed in this work for cases with larger amount of absorbent. The general comparison between the different methods shows * Note that a composition error of the B 4 C was introduced in the numerical results of [5] and [3]. This error was intentionally reproduced here to be able to compare our results with those found in the literature. * RMC used as reference: k ∞ and ∆k ∞ = (k HEBE ∞ − k ∞ ) * 10 5 underneath that they can be in agreement as long as the amount of absorbents (either fuel or poison) is not too large, and that the equivalent fraction is chosen adequately for the case. This observation represents actually a drawback to the theory that a constant fraction should be taken. More work is in progress to identify the source of the discrepancy. However, since it increases with the fuel or burnable absorber content, we can suppose that the approximation of independent particles may not be valid in those cases.