D continuum percolation approach and its application to lava-like fuel-containing materials behaviour forecast

The paper is devoted to the theoretical study of elementary permeable objects percolation and its application to real physical objects. Spheres and isotropic oriented capped sticks were chosen as elementary geometrical objects for percolation simulation, physically adequate for radiation defects behaviour description in brittle dielectrics, particularly in the so-called Lava-like Fuel Containing Materials (LFCM), where it effects their mechanical steadiness. LFCMs is high-radioactive glass, which was formed during active stage of well-known heavy nuclear accident, that occurred at Chornobyl nuclear facility in 1986. Physical processes taking place in the materials are of great practical interest. Furthermore, when applying percolation models to LFCM objects, an approximate behaviour forecast can be created. From the results of simulation, it appears that physical properties of the LFCM should drastically change within in the period of 2015÷2045 calendar years, depending on variations in nuclear fuel content.


Introduction
The percolation approach is widely used in solid-state physics, particularly for its description as a unified system, containing inclusions of another structure.To these may be referred composite materials, glass-ceramics, and complex structures with radiation damages in its volume, and many others.Roughly, a generalized three-dimensional percolation problem can be formulated in the following manner: there are some random non-homogeneous media, for example, stone with randomly distributed pores, and there is some liquid (water) or gas from one side of the stone, which may percolate inside and through the pores.And there arises a question, concerning the conditions (concentration, size, form, positions of the pores) under which the substance will propagate (percolate) from one side to the other.Such conditions are called as critical conditions or critical points.Those sorts of tasks are known as percolation problems.Many static and dynamic physical properties of multiphase materials and systems are determined by microscopic and macroscopic spatial distribution of their phases.When the system or material is near a critical point or starts to undergo a structural phase transition, the properties of such a system are best modelled by continuum percolation with objects of a corresponding shape and size and are properly (in most cases randomly) distributed, in contrast to site or bond percolation where sites or bonds are in a discrete lattice and randomly occupied.In most systems inclusions of another structure are physically small and randomly distributed, which makes percolation the best tool for modelling the phase transitions and critical behaviour in unified systems, containing the inclusions.
Percolation problems of this kind are known as three-dimensional (3D) continuum percolation of hard-core or soft-core (permeable) geometrical objects and they were the object of persistent research in the 80's [1][2][3][4][5][6].Due to computational complexity such percolation tasks, however, are much less known than the well-known lattice percolation tasks, though the progress in computational techniques in the last decade has made it possible to achieve significant results in solving some of them [2][3][4]7,8].
The percolation approach is widely applied to condensed matter problems solved in the way of simulation.Common examples are doped materials [9] and polymers [10].The simulation of etching process in natural application for percolation simulations as well [11].The present work is devoted to the solution of some special percolation problems in order to characterize critical behaviour resulting from volume-generated radiation damages in condensed matter.In particular, percolation models of permeable spheres and sticks are used; the real physical object for application is LFCM (Lava-like Fuel-Containing Materials), relating to the class of brittle materials; radiation defects accumulation, as it will be shown below, can lead to critical change in its mechanical properties and behaviour in future.From the practical point of view, LFCM mechanical steadiness to external impacts is of especial interest, because LFCM mechanical destruction and transformation into high-dispersive state, regarding its high specific radioactivity, is a potentially dangerous event.

Percolation problem on soft spheres
The mathematical formulation of the continuum percolation problem on the soft spheres is as follows: there is a cube; it is filled with small (relative to the cube size) permeable spheres, their centres (i.e.nodes) being randomly distributed.The critical concentration is to be reached when it is possible to find a way through spheres from one side of the cube to the opposite side (so percolation is reached).Such a task is known as three-dimensional percolation problem on random sites.It is obvious that percolation is attained when infinite (in terms of percolation theory) cluster from the sites is formed.The most suitable parameter to characterize the percolation process independently from the linear dimensions, concentration and even the form of objects on which the percolation is realized is average bond connections per node at percolation threshold B c : Here N is the number of nodes, V is volume of the simulation cube, which is taken unitary for simplicity, v ex is an excluded geometric figure volume.The excluded volume is defined as a volume surrounding and including a given geometric percolation figure (sphere in the present case), which is inaccessible to another figure [12].For spheres, this is a sphere of double size, where it is easy to get the B c for soft spheres: where r c is the critical spheres radius at which percolation is realised.The number of nodes should be fairly large due to a systematic error (so-called finite-size effects) increasing with the size of spheres, which slowly decreases with an increase of the number of nodes.To obtain the threshold B c in a "direct" way, finite-size effects should be excluded, which requires zero percolation radius, and, consequently, an infinite number of percolation nodes.Nonetheless, the exact percolation threshold with uncertainty as small as is required could be obtained by excluding a systematic error of finitesize effect in the following way: a dependence of percolation threshold versus spheres radius should be built, then approximated and extrapolated with power function.Extrapolation of B c value to zero percolation radius and a corresponding infinite number of nodes can be considered as a "true" value for three-dimensional percolation threshold.The reverse of power in approximated power function is commonly called a critical index ν [13].
A lot of algorithms for solving percolation problems have been developed [3,4,7,8,12,13].On having revised them, the author found out that they lacked speed and were doing a lot of unnecessary calculations in the process, repeating them over and over again, and even the fastest algorithms [4] are not free from certain drawbacks.Considering that the author has developed a special algorithm for percolation on random sites, and for the first time it was tested for percolation on spheres.Herein below there is a brief description of the main steps of the algorithm.
1. Random distribution of the nodes within the simulated cubic region is made, the nodes positions are set to the proper database with X,Y and Z axes oriented along the cube axes; 2. The positions of the spheres are sorted in the positions database index in order to quickly figure out the neighbouring spheres around the chosen one.For example, in one index they are sorted by their X coordinate in ascending manner.It means that at the stage of defining distances between the nodes one has to take into account the nodes only within a known range of corresponding index.A similar procedure of indexing is done for Y and Z coordinates.
3. For each node a distance to nearby nodes is determined.
Step 2 makes the process much faster.For each separate node, its neighbours within a certain range on X, Y, and Z axes can be quickly determined.All neighbours are automatically placed together in the database, because they are sorted by corresponding X,Y or Z coordinates.
Let us explain the item 3 more in detail: 3.1 First, the program gets coordinates of the current node and its position in the database.
Next, in the database, the first index array of all the nodes is sorted by X coordinate with the current node as starting point.Then, the program goes up on the sorted array starting from the current node position, picking the numbers of the nodes that are close to the current one by X coordinate.This process stops when the difference between X coordinate of the current node and the other node becomes too large, i. e., more than some limiting distance (more than the distance at which percolation surely occurs).The same procedure is performed, but the array is run down by X coordinate.So, in the end we got recorded an array of nodes which are close to the current node by X coordinate.The same procedure is performed for Y and Z coordinates, and only those nodes are checked for distance that are close enough to X, Y, and Z coordinates simultaneously.By choosing the correct limiting radius, the number of nodes finally checked for the distance can be reduced to the average of 4-5, instead of running the full array of nodes required to be checked for distance.Therefore, time consuming radius and distance calculations are replaced by a rather quick picking up the node numbers from the database.It runs much faster, and, in addition, its calculation time is roughly proportional to N 4/3 , or even to N for heavilyoptimized case, with time for initial sorting and indexing roughly proportional to N lnN (quick sorting); N here is the number of nodes.For comparison, calculation time is proportional to N 2 if one uses a "classical" approach without sorting.Therefore, for huge node numbers ( > 100, 000) the calculation scheme proposed can be very effective.Actually, for N > 500, 000, the time required for calculation of relative positions was less than all other steps of algorithm taken together.Thus, for a large ( > 10, 000) number of nodes, the preliminary indexing makes sense despite the time required for the indexing procedure.All the procedures at the start take a lot of memory for indexing, for arrays, etc., but for the current situation, memory was sacrificed for the sake of speed.

3.2
The distances are calculated and recorded for corresponding pairs of nodes.After the procedure, the database of indexes and even positions can be eliminated: the operation turns from three-dimensional problem into simple logical pathfinding problem.
4. After turning the assembly of nodes from spatial structure into a connected graph, a special pathfinding algorythm, which can be classified as conditional propagation pathfinding algorithm with some author's modification, is applied.It finds a path through the nodes from one edge of the cube to the opposite one with optimization by a minimum distance between the nodes on the path.Technically, this is realized in the following way: every node has a specific parameter assigned: an individual "percolation radius", e.g., characteristic radius of spheres, at which the current node will be percolated.
At the starting iteration, all the nodes have the "percolation radius" equal to their distance from the starting edge of the cube, e.g., they percolate when the sphere centred on the node "touches" the starting wall of the cube.Then, each pair of nodes is checked in the following way: if the maximum of values of "percolation radius" of the first node and half the distance between nodes is less than "percolation radius" of the second node, the percolation radius of the second node becomes equal to the maximum of the two values.This process is illustrated schematically in figure 1.The left-hand node in figure 1 has "percolation radius" which is equal, naturally, to the distance from the node to the cube wall, which is a starting one for pathfinding.The node will percolate if the sphere radius is sufficient to touch the starting wall.The left and middle nodes will touch each other if the percolation radius reaches half the distance between the nodes, and "percolation radius" for the left node is less than half the distance between the left and the middle nodes.So, the middle node will percolate at the sphere radius equal to half the distance between the left and the middle node.Thus, "percolation radius" of the middle node is assigned as being equal to half the distance between the left and the middle node.The right node can percolate from the middle node only, and its "percolation radius" is equal to the radius of the middle node, because half-distance between the nodes is less than "percolation radius" of the middle node.So, the right node will percolate at the sphere radius equal to the radius of the middle node.
In the process of successive steps, if the maximum value obtained is larger than "percolation radius" of the current node, "percolation radius" naturally does not get assigned to that maximum.
Therefore, in successive procedure of "percolation radius" re-assignment, the modified pathfinding algorithm, in a few iterations, re-assigns the "percolation radius" to every node.As the final result of iterations, all "percolation radius" values became minimal, e.g., if the sphere radius is less than "percolation radius" the node will certainly will not be percolated.
Thus, to obtain the radius at percolation threshold one should get the "percolation radius" from the nodes at the remote side of the cube only.Then, knowing N and the resulting percolation radius r c , one can easily calculate B c .
The algorithm is rather general and can be applied to a system of virtually any number of dimensions and geometric shapes.This algorithm is also very convenient in defining the dependence of infinite cluster power on the sphere radius.For accounting one should only build a dependence on the numbers of spheres whose "percolation radius" is equal or less than the sphere radius.
The method turned out to be rather robust and reliable.So, the percolation threshold determination procedure after initial tests underwent slight changes: several initial percolating nodes were placed in the centre of the cube with "percolation radius" property assigned to zero, and, after a few iterations, it is easy to get a dependence of cluster power on the radius and thus to obtain the percolation threshold by power-law approximation.
For each given number of nodes, a few thousand variants of node positions were taken initially for every simulation, B c calculated for every variant, and, in the end, they were averaged.The spread of percolation thresholds for different realizations can be estimated from 2: the probability of percolation realization is essentially an integral of percolation event probability density.Subsequently, a steeper curve corresponds to a narrower distribution.
As a conclusion, this algorithm has the following advantages: no calculation time wasting on adding another node and on re-calculating all the node assembly after every increment; only nearby nodes are taken into account which are simply taken from the database without any additional calculations; r c and B c is found in just one pathfinding session.Compared to common algorithms, [1,2,4], in practice this is about 10× increase in productivity for 100,000 nodes, further increasing with the nodes number.Con is rather obvious: the algorithm is correctly applicable to fully permeable objects only.
For objects of non-spherical shape, one should understand property of "percolation radius" as a scaling factor applied to all geometrical figures, provided their position kept constant and figures are convex.In such a fashion, one can easily apply this algorithm even to the objects with various but proportional sizes, and even to objects of different shape.
The dependence of probability of percolation (ratio of events where percolation occurred to the total number of events) on B (the average number of bond connections per node which is proportional to density) is presented in figure 2. The curves were calculated for different node numbers to evaluate the error given by finite concentration.For every N , finite size effect suggests different B c as well as r c .In figure 3 the plot of B c versus r c is presented, each point representing the results for its own N .Error bars are actually about the size of a dot and thus they are not shown.The dependence very well fits to the power type.One can extrapolate the true B c at the point of r c = 0, and at this point B c is 2.734 ± 0.005.The power is 1.11 ± 0.02.It corresponds to the percolation critical index of ν = 0.89 ± 0.02, similar to the results established in [13].Compared to results established earlier, the current realization presents a more efficient algorithm only, which, using appropriate calculational power, can be used for a more precise determination of percolation parameters.

Percolation task on soft sticks
The mathematical formulation of the task is as follows: there is a cube; it is randomly and isotropically filled with fully permeable objects, having the form of capped cylinders.The percolation is realized on the sticks from one wall of the cube to the other.The first solutions of such a task are in [1,12].The solution revealed the analytical approach to the solution of such a task, known as an excluded volume rule.However, their simulation was not quite correct, as it is revealed in further work [2] due to the flaw in the mathematical model built: the distribution of the stick's direction was not truly isotropic.
We have excluded the flaw in our model in principle by generating direction vectors in the following way: a random radius vector was generated within the volume-centred cube with edge length twice the unitary.If the length of the radius vector exceeded the unitary length, thus taking place out of sphere with unitary radius, the vector was generating again until a vector of length no more than unitary was obtained.So, the endpoints of the vectors distributed uniformly inside the sphere with the starting point in its centre and, therefore, isotropic.Afterwards, the vector was divided on its length.As a final result, vectors of unitary length and isotropic distribution were obtained.The process is roughly illustrated in figure 4a: the radius-vector shorter than the unitary, stays, longer than the unitary rejected.All that is left to do is to normalize them to have a the unitary length.
(b) If one takes all the radius-vectors within the cube, for example, the ones directed to the corners of the cube it is more probable that the distribution will not be isotropic.In work [2] the distribution was on the spherical coordinates, so the random distribution looked like in figure 4b i. e., denser at the poles of spherical coordinates.The method described above eliminates the flaws of that kind in principle.The computational algorithm is similar to the one used with permeable spheres; the radius of cylinders, at which pair of capped cylinders at a given position and directions came in contact while maintaining R/L ratio, was taken as "percolation radius" property.Here L is the cylinder length and R is the cylinder radius, as shown in figure 5. Apart from adding the property of direction and somewhat complex algorithm of determining "percolation radius" to every pair of nodes, percolation algorithm required no changes in comparison with percolation on spheres.The simulations were run with different ratio of the capped cylinders radius R to their length L. For each R/L value, the simulation of percolation was performed with different concentration of the capped sticks N .Then, like in our percolation with spheres, a power extrapolation was made as shown in figure 7 and the dependence of extrapolation on R/L ratio is presented in figure 6.Thus, each point on the graph is an effective approximation of percolation threshold by different numbers of objects in volume.
In the work we have made a slight advance in R/L ratio in comparison with [2], so we can observe that power-law at low R/L maintains up to R/L equal to 0.001 and the power is equal to 0.58 ± 0.02.In general, the curve presented in figure 6 lies slightly below the curve presented in [2], which can be explained by more accurate determination of B c .

Application to the real material: LFCM behaviour forecast
LFCM (Lava-like Fuel-Containing Materials) are the materials formed as a result of the wellknown heavy nuclear accident in Chornobyl in 1986.As it is accepted for now [14], the LFCM are glass ceramic alkaline-earth silicate compositions (devitrified glasses) containing 5-10 mass percent of irradiated uranium nuclear fuel in their volume in company with high-radioactive fission products.The review of the main LFCM physical properties is presented in [14].Approximately 1,000 tons of this material, accompanied both by the other core debris and the building constructions of the destroyed Chornobyl NPP 4-th unit, is located in the so-called "Shelter" site, which had been erected quickly after 1986 accident as a forced measure directed to primary prevention of radioactive substances dissemination in the environment.The "Shelter" site is not a hermetically sealed construction and there are no doubts that a large quantity of high-radioactive compositions, such as LFCM and irradiated nuclear fuel itself do have a direct air contact with the environment.LFCM themselves look like the coloured glass ceramics and one can classify them like brown (8-10% fuel content), multicoloured (6-7%) and black (4% fuel as usually).These materials do have a lot of interesting properties, among which there is a spontaneous self-sputtering of its surface [15].The percolation models presented above can be successfully applied to the ensemble of radiation damages behaviour modelling.The LFCM do have disordered regions (DR) of spherical form in its volume, induced by heavy recoil nuclei [16].The DR spatial distribution is quite random, and the neighbouring DR can be overlapped in a great measure, which just coincides with the model of permeable spheres.Therefore, the model of percolating permeable spheres looks physically adequate.If we apply the above results of the modelling to LFCM, we will acquire the results given below.First of all, we can establish the approximate time when the infinite cluster of disordered regions will be formed.Taking into account the variations in LFCM specific radioactivity and rate of defect formation, one can approximately determine the DR concentration and, therefore, the calendar year, when the infinite cluster will be formed.Taking the DR diameter to be 25 nm [16], knowing formula (2) and B c for spheres, calculated earlier, we can derive the critical concentration of DR: (2.5 ± 0.5) • 10 16 cm −3 .The most of the error arises from variations in DR diameter, which, as it was established in [16], is about 25 ± 2nm.Considering the rate of defects generation (8 ± 2) • 10 14 cm −3 /year for 10% fuel content [16], the formation of such a concentration to the calendar year of 2015 ÷ 2045 depending the variations in alpha-active nuclides concentration.Probability of such a percolation event in the calendar year for LFCM containing 10% of fuel is presented in figure 8.
Finite size effects expressed in different curves for different numbers of defects in this case do have its own physical sense: when there is a huge block of LFCM, a curve for the greatest number of defects should be taken into account.With time, however, LFCM will undergo further destruction into individual dust particles of a smaller grade, containing lesser number of defects, so the curve for 10,000 or even for 1,000 should be appropriate.Thus, in accordance with figure 8, the process is likely to develop in the following way: initially, large (macroscopic) LFCM fragment will be fractured into cluster particles of the order of tens micron grade.Then (this will occur fast enough) such fragments will fall apart into "elementary" particles with the grade of one or a few DR.This final stage of destruction means conversion of the whole LFCM volume into submicronic dust, just similar to that observed in experiments on LFCM surface self-sputtering [15].
The above description of fracture process is based on fundamental concepts of fracture mechanics [17].The volume fracture of brittle materials starts with cracks development, which propagate along the border, splitting the areas with different mechanical properties.Such a kind of border one may identify as a cluster edge.This version has a confirmation in experiments, when spontaneous dust productivity of brittle dielectrics (LFCM) was observed: the shape of dust particles was identified as clustered disordered regions, i. e., mechanical destruction of the surface occurred along the borders dividing DR from the surrounding material [15].
With the increase of the density of defects, the mean cluster spanning size increases as well, obviously leading to an increase in the average lengths of the cracks.At this stage, however, there will be no radical changes in material properties, since the crack propagation will be limited by the nearest disordered region.The situation will change drastically, however, when the infinite cluster is formed.
In light of the fact that cracks appear at the borders of the disordered regions, there is no doubt that the infinite cluster of disordered regions formation in LFCM means the formation of the main crack along the infinite cluster border, penetrating the macroscopic fragment of material thorough as a whole.Considering percolation properties of the infinite cluster (spanning all along the system, occupying most of the volume), LFCM after percolation threshold will likely no longer possess any macroscopic mechanical properties: it will be reduced to fine grade sand.The approximate LFCM volume dust productivity increase one can roughly estimate as 20× minimum.The further additional changes in the material structure may occur when the area of the clustered DR surface reaches its maximum -this will cause the leavings of LFCM to be especially susceptible to mechanical loads and other external factors.The clustered DR surface parameter is of great interest due to its effect on mechanical steadiness of materials beyond the critical points.Therefore, due to the practical importance of the mentioned parameter, the dependence of such an area per unit volume was calculated; corresponding results are shown in the figure 9.The task was calculated accepting DR diameter to be 25 nm, like it takes place in LFCM.One can see from figure 9 that specific border surface will reach its maximum at 6.5 • 10 16 cm −3 DR concentration, which corresponds to the calendar years of 2050 ÷ 2080.
If one tries to apply the theory of Coulomb explosion [18] and thermal spikes [19] to the LFCM, we will get the percolation task on permeable sticks, since Coulomb explosion forms stick-like heavily disordered regions in LFCM.Considering density of sticks as equal to the DR density, application of the percolation based on the soft sticks calculation results to the example yields the percolation threshold at the concentration of 9.4•10 19 cm −3 , which is much larger than the numbers for the model of permeable spheres, but much less than the known numbers of steadiness for silicate glass under alpha-irradiation.Such a concentration corresponds to 110,000 years of internal irradiation action.Nonetheless, such a structure of intersecting percolation stick can effectively serve as mechanical stress concentrators, thus making the material investigated rather susceptible to mechanical loads.

Conclusions
The 3D percolation problem on random sites with permeable spheres or sticks as percolating objects is an appropriate way of stimulating the behaviour of radiation damages ensemble in LFCM.
Physical considerations, coming from the basics of the fracture mechanics of brittle materials, make it possible to conclude that LFCM will eventually mechanically degrade into highly dispersed state.The current results of computer simulations based on percolation approach indicate that the potentially dangerous process of LFCM destruction will be urgent in observable future.
Analysis of the results obtained indicates that the destruction of LFCM will undergo in two stages: first, LFCM will be transformed into fine grain sand (having typical grade of a few microns); second, the sand will quickly turn into submicronic dust.Quantitative calculation results allow to predict that for brown LFCM with most of fuel content(8%) the destruction will be actual starting 2015-2020; for black LFCM with less fuel content(4%) this may be actual after year 2045.

Figure 2 .
Figure 2. Dependence of percolation probability W on the average bond connections.Calculations for different node concentrations per simulation box were performed from 1,000 up to 100,000 spheres.

Figure 3 .
Figure 3. Plot of percolation threshold Bc versus critical radius of spheres rc.The simulation cube is of unitary size.

Figure 4 .
Figure 4. Illustration for isotropic vector generation procedure.On the left you can see (a) the choice of isotropic radius-vector and on the right (b) the generated radius-vectors with random spherical coordinates.

Figure 6 .
Figure 6.S = Bc − 1 as a function of the R/L aspect ratios of sticks presented on the left graph.On the right graph, the dependence of S on R/L for the range of R/L < 0.01 is shown in the log-log scale.The straight line is the power function approximation shows the power of 0.58 ± 0.02.

Figure 7 .
Figure 7. Dependence of percolation threshold on sticks grade.Triangles correspond to R/L = 0.02, circles correspond to R/L = 0.01 and squares correspond to R/L = 0.005 Percolation is performed inside a cube of unitary size.

Figure 8 .
Figure 8. Probability of percolation event versus calendar year plot for LFCM containing 10% of fuel.

Figure 9 .
Figure 9. Plot of total spheres clusters specific surface versus defects density.