PRELIMINARY INVESTIGATIONS OF TRANSPORT IN HETEROGENEOUS RANDOM MEDIA

Random media emerge in several applications in reactor physics and safety analysis. Most often, models of stochastic media assume spatial homogeneity, whereas real-world complex materials, such as fuel chunks resulting from core degradation, typically display apparent heterogeneities. In a series of previous works, we have shown that stochastic tessellations can be successfully used in order to describe the material properties of several classes of random media. In this paper we extend these results to the case of heterogeneous random media by using Voronoi tessellations with space-dependent seed distributions, allowing for spatial gradients.


INTRODUCTION
Stochastic media are key to many applications in reactor physics [1], encompassing the evaluation of the re-criticality probability following severe accident with fuel degradation [2,3], the analysis of neutron absorbers with dispersed poison grains [4] or MOX fuels with Pu-rich agglomerates [5]. In this context, one is typically interested in determining the angular particle flux ϕ , averaged over the ensemble of possible material realizations corresponding to the underlying disorder model. We will consider the case of binary stochastic mixtures, where at each spatial location two immiscible materials (say α or β) are assumed to be randomly present [1]. Reference solutions for the ensemble-averaged neutron flux ϕ can be obtained by using the quenched-disorder approach [1]. A set of material realizations is first sampled from a given statistical model; then, the corresponding transport problem is solved for each realization (assumed to be frozen), and the quantities of interest are computed; finally, the ensemble averages or the full distribution are determined with respect to the ensemble [1]. This procedure is exact, but highly time-consuming, and can thus be used in order to validate approximate but faster models based on the annealed disorder approach: a single transport calculation with effective particle displacement laws is performed in a 'homogenized' medium, at the expense of neglecting the contributions of disorder-induced spatial correlations [1,6,7]. For both approaches, the properties of neutron transport within random media are intimately related to the average chord length Λ i of the underlying stochastic geometries, i.e., the typical length of the intersection of a given line with the random interfaces of the material chunks of kind i = α, β [1]. When the neutron mean free path 1/Σ t,i for each material i is much larger than Λ i , each particle flight will experience a large number of material chunks, and we fall thus within the atomic mix regime: we can safely assume that the medium is composed of a homogeneous material, whose properties are obtained by averaging the properties of each random material, weighted by their respective probability p i of occurrence. In the opposite case, there is a non-trivial interplay between the flight lengths and the material chunk size.
In a series of recent papers, we have shown that stochastic tessellations, which are probabilistic models where a given spatial region is randomly partitioned into convex cells [8,9], can be successfully used in order to generate the random material realizations for the quenched-disorder approach [3,[10][11][12]. We have examined several classes of tessellations, including Poisson and Poisson-Voronoi, and we have computed reference solutions for non-multiplying and multiplying benchmark configurations [3,11]. In our previous investigations, we have assumed that statistical properties of the tessellations (and in particular the associated average chord lengths) were the same at any spatial position. In several real-world applications, one might encounter random media displaying a gradient along a given direction, typically close to the boundaries [13,14]. Such situations are common, e.g., in material science, where layers with structural inhomogeneities can either result from manufacturing or from unwanted stratification due to gravity or pressure [14]. Non-stationary stochastic tessellations have received so far only a limited attention [14,15]. Recent works have addressed particle transport in one-dimensional spatially heterogeneous annealeddisorder models [16]. In this paper we will present a preliminary investigation of neutron propagation in three-dimensional stochastic tessellations with a spatial gradient. For this purpose, we will suitably modify a Poisson-Voronoi tessellation model by introducing position-dependent random 'seeds' * , which can be explicitly sampled by Monte Carlo methods. Then, we will probe the properties of neutron transport within these heterogeneous random media for some simple nonmultiplying benchmark configurations and we will compare our numerical findings to the case of spatially homogeneous stochastic tessellations.

HETEROGENEOUS POISSON-VORONOI TESSELLATIONS
Voronoi tessellations are a prototype process for isotropic random division of space, often in connection with material fragmentation models [8]. A portion of a space is decomposed into polyhedral cells by a partitioning process based on a set of random points, called 'seeds'. The Voronoi decomposition is obtained from this set of seeds by applying the following deterministic procedure: each seed is associated with a (convex) Voronoi cell, defined as the set of points which are nearer to this seed than to any other seed. For spatially homogeneous Poisson-Voronoi tessellations, the seeds obey a homogeneous Poisson point process with constant density τ [9]. For the sake of simplicity, we will assume that the mean variability of the stochastic medium is along the x direction, whereas the probability distribution of the spatial realizations is invariant by translations parallel to the (y, z) plane (for an illustration, see Fig. 1).
Gradient structures can be introduced in Poisson-Voronoi tessellations by sampling the seeds from a non-homogeneous Poisson process with intensity τ (x), over 0 < x < L x , L x being the lin- For each law, we show the plane section Oxy of the tessellation at z = 0 before the coloring procedure (left), the three-dimensional visualization before the coloring procedure (middle) and the three-dimensional visualization after the coloring procedure with probability p α = 0.1 of assigning a red label (right).
ear size of the random medium with respect to the x direction. By building upon the computer code that we have recently implemented for stationary Poisson-Voronoi tessellations [17], general laws τ (x) can be taken into account with minor modifications. The first step is the sampling of the number of seeds to be introduced in the tessellation from a Poisson distribution with mean τ w = τ (x, y, z)dxdydz. In the special case of gradient structures along x-axis within a box For each seed, a trial position (x, y, z) is uniformly chosen within the box: the position is accepted with probability τ (x, y, z)/τ w ; otherwise, the position is rejected and has to be re-sampled until acceptance.

STATISTICAL PROPERTIES OF HETEROGENEOUS POISSON-VORONOI TES-SELLATIONS
The average cell volume V and the average chord length Λ of a tessellation deeply affect transport properties † . For three-dimensional stationary Poisson-Voronoi tessellations with constant τ , their † Here Λ is the average chord length of the tessellation, prior to assigning the material properties. For binary mixtures, Λ i for the colored tessellation is then obtained as Λ i = Λ/(1 − p i ).
analytical expressions are known [9], and read which clearly depend on the tessellation density τ . For gradient structures with position-dependent τ (x), in a few cases analytical formulas can be derived, but they are generally expressed as integrals that can be hardly exploited numerically [14]. Nonetheless, it has been shown that V (x) and Λ(x) can be approximated to a high degree of accuracy by remarkably simple expressions [14], namely, We have numerically verified the formulas given in Eq.
and finally sine-like variation, with h) Λ(x) = 0.1 + 0.00064x 2 (x − 10) 2 and i) Λ(x) = 0.5 − 0.00064x 2 (x − 10) 2 . Some realizations are illustrated in Fig. 1. For each law τ (x) we have generated a collection of realizations and we have then estimated the ensemble averages V (x) and Λ(x) by probing the volumes of the cells and the length of random chords crossing the tessellations at each position x. The results corresponding to these gradient laws τ (x) are shown in Fig. 2 for the computed average chord length (left) and the average volume (right), displayed as a function of x. For each observable, simulation results have been obtained by two ways: the first option consists in sampling the seeds between x min = 0 and x max = L x ; the second option consists in extending the sampling of the seeds along the x-axis by taking x min = −1 and x max = L x + 1, in order to avoid finite-size effects at the boundaries of the box. When geometries are produced with the first method, strong finite-size effects occur in the neighbourhood of the boundaries (x = 0 and x = L x ). The impact of finite-size effects is greatly reduced by using the second technique, as shown in Fig. 2, especially for the chord lengths. This is why, in the following, we will sample seeds in the extended domain and extract the portion of the geometry between x = 0 and x = L x for statistical analysis on for the transport problems. Generally speaking, the agreement between the Monte Carlo simulation findings and the approximation provided by the analytical formulas in Eq. (2) is fairly good for the volume and for the chord length. As expected, the discontinuous cases (f and g) yield the worst agreement in the neighbourhood of the step variation. The discrepancy with respect to the analytical curves stem from both the formulas in Eq. (2) being an approximation (albeit accurate) and from finite-size effects. As a general trend, the agreement decreases with increasing Λ(x). In particular, for the case corresponding to the largest large chord length, the effective volume is larger than predicted by Eq. (2), in spite of the presence of truncated polyhedra (which should reduce the average volume); this behaviour was previously remarked in [18].

PARTICLE TRANSPORT IN BENCHMARK CONFIGURATIONS
In order to probe the behaviour of neutron transport within gradient structures, we have selected the following configuration: particles undergo single-speed transport with scattering and absorption in a box of side L x = 10, where reflecting boundary conditions are applied on all faces, except the two orthogonal to the x axis and located at x = 0 and x = L x . A uniform and isotropic source is applied on the face at x = 0. The stochastic medium is a binary mixture composed of an absorbing material, denoted α, and a scattering material, denoted β; this mixture obeys Poisson-Voronoi statistics, with density τ (x) obeying one of the laws described above; the coloring probability of the mixture is p α = 0.1. We consider two sets of material cross sections (in arbitrary units): either Σ α t = Σ α a = 10 and Σ β t = Σ β s = 1; or material α is unchanged and material β has Σ β t = Σ β s = 10 (this latter case is chosen in order to assess the interplay of the cross sections and the chord lengths). The number of replicas for ensemble averages varies between 200 and 1000, depending on the law and on the material compositions. For each replica, particle transport is simulated by using the Monte Carlo code TRIPOLI-4 R , developed at CEA.
Simulation results are shown in Fig. 3. For both sets of material cross sections, the impact of the laws τ (x) on the particle flux ϕ(x) is clearly visible. The limit behaviour is naturally imposed by the extreme case of constant minimum Λ = 0.1 and constant maximum Λ = 0.5. In all tested configurations, the flux ϕ(x) associated to gradient structures cannot be reasonably well approximated by using the a constant value of Λ equal to the average between the limit cases (law b, with Λ = 0.3). Due to the particular choices of τ (x) given above, close to the source position the flux corresponding to heterogeneous geometries initially follows the one corresponding to the limit cases. Deviations increase as τ (x) progressively deviates from the value at x = 0: this is particularly apparent for the case of a step heterogeneity at the center of the box (Fig. 3, center). For the case where Σ β s = 1, the pairs of laws τ (x) that are symmetric with respect to each other (d and e, f and g) seem to lead to very similar flux values at the end of the box (x = L x ); this might be the case also for the case where Σ β s = 10, but the flux is decreasing rapidly with x and statistical uncertainties prevent from drawing solid conclusions. For reference, the atomic mix regime (which would correspond to the limit of very small Λ i for both materials) is also shown in Fig. 3.

CONCLUSIONS
In this work we have examined a preliminary investigation of particle transport through threedimensional random media with spatial gradients. For this purpose, we have used a quenched disorder approach based on Poisson-Voronoi tessellations. The statistical features of these geometries and the properties of the ensemble-averaged particle flux have been examined for a few relevant gradient profiles within a box. It would be interesting to compare the reference results obtained in this work to annealed disorder models. Future investigations will also address the case of heterogeneous Poisson hyperplane tessellations, which form a prominent class of random media and are commonly used in several applications.  ϕ(x) x Figure 3: Ensemble-averaged spatial flux ϕ(x) , in a box of size L x = 10, for different densities τ (x). Green inverted triangles correspond to law a, black diamonds to law b and orange triangles to law c. Red circles correspond to law d (top), law f (center) and to law h (bottom). Blue circles correspond to law e (top), law g (center) and law i (bottom). Left. Σ β s = 1. Right. Σ β s = 10. For comparison, grey stars represent the atomic mix case.