Stochastic generation of virtual air pores in granular materials

A computational method is described for the generation of virtual air pores with randomized features in granular materials. The method is based on the creation of a stack of two dimensional stochastically generated domains of packed virtual aggregate particles that are converted to three dimensions and made to intersected with one another. The three dimensional structure that is created is then sampled with an algorithm that detects the void space left between the intersected particles, which corresponds to the air void volume in real materials. This allows the generation of a map of the previously generated three dimensional model that can be used to analyse the topology of the void channels. The isotropy of the samples is here discussed and analysed. The air void size distribution in all the virtual samples generated in this study is described with the Weibull distribution and the goodness of fit is successfully evaluated with the Kolmogorov–Smirnov test. The specific surface of the virtual samples is also successfully compared to that of real samples. The results show that a stochastic approach to the generation of virtual granular materials based only on geometric principles is feasible and provides realistic results.


Introduction
The study of porosity in granular materials is important to understand the behaviour of a number of physical properties, e.g.Young's modulus [23], strength [23], noisiness [28], and infiltration rate [10].In addition to these, in the field of civil engineering, porosity is related to other properties, such as the durability of pavements [20] or their resistance to rutting, fatigue cracking, and low temperature cracking [17].Therefore, it is necessary to take porosity into account when designing granular materials for various applications.
The study of porosity, however, is usually pursued by performing X-ray CT scans of the specimens under analysis [12], which is an expensive and time consuming process.Due to their cost, X-ray CT scanners are not accessible to all, thus, computational methods were developed to reproduce the internal structure of materials.The models available in the current literature take a variety of approaches and range from the use of the discrete elements method (DEM) to analyse the particles in a granular material [4,9] to the use of the percolation method to describe the void space [11,24].Moreover, in [5] the authors show that three dimensional particulate models obtained with DEM can be further elaborated to perform an analysis of the void space, too.The current models describing porous media are all based on advanced theoretical principles, therefore, their implementation is complex, unless commercial software is used [21] or the researchers are expert in mathematical modelling.
In [6] a new model based on geometry only is described and analysed by the present authors.The method was developed in order to allow a theoretically-based and simple way to generate granular material virtual representations.In [6], the authors proved that the method provides realistic results, however, a further investigation on the features of the void space is necessary, mainly because the 3D model generated previously only provided a proof of concept and was obtained through a manual trial and error method.
In this paper, we describe the automation of the mathematical method presented in [6] and analyse the void patterns that it generates, providing topological information on the virtual samples that are created in terms of isotropy and coordination number.In addition, the void size distribution of the virtual porous media is analysed.
The main aim of this paper is providing a first validation of this newly developed model for the generation of air pores in granular materials by comparing the virtual air voids to the air voids found in real samples.

Generation of 3D virtual granular materials 2.1 Packing of virtual granular elements in a 2D domain
In [6] a method to generate virtual asphalt samples is described.The method starts with the generation of 2D planes of virtual particles based on the biological mechanism of contact inhibition and can be addressed as a packing algorithm.A number of virtual particles (circles or ellipses) are seeded in random points in an arbitrary subset of R 2 and grown until they reach a fixed maximum radius, reach the border of the domain, or touch another particle.Therefore, if we assume that at any point of the analysis all particles respect the just mentioned criteria and they all are circles, the conditions for the growth are: where r i is the radius of the ith particle, Δr is the fixed radius increase, r max is the maximum permitted particle radius, C r i +Δr is the set of points defining each circumference after the radius increase, A ⊂ R 2 is the arbitrary domain of interest, and d is the Euclidean distance (2D) between the centre of particle i from the centres of all other particles.The role of the domain A is further discussed in Sect.5.1.A set of conditions similar to relation (1) can be developed in the case of ellipses, too.For further details on the packing method for ellipses, see [6].
The randomization of the coordinates of the centres of the particles in the selected portion of the x y plane [set A in relation (1)] is performed with a random number generator using the standard uniform distribution [6].This distribution is described by the following probability density function: Fig. 1 Sample domain generated with the algorithm described in [6] where a = 0 and b = 1.It is important to mention that the centres of growth in our method form a simple binomial point process in the compact set M = [a, b] 2 [7].The property of simplicity refers to the fact that with probability 1 no points of the process may coincide [3].In addition, since the shape chosen for the generation of 2D sections in this study is a circle in the plane x y, rejection sampling is applied [7]: a sequence of n independent uniform random points is generated in M until a point falls in the set W ⊂ M and defined by the equation where c x = c y = b/2 and r = 0.5.This operation is then repeated until a satisfactory (user-defined) number of centres of growth is reached.Furthermore, scale changes are applied to the centres in order to reach the desired diameter for the virtual cross sections.The shape of the subset W was set as a circle (Fig. 1) in order to allow a comparison with X-ray CT scans of commonly used asphalt samples.In the case of a rectangular domain, the set W would be defined as Even if the centres are seeded according to a uniform distribution, the radii of the particles composing the domains generated with our method do not show the same behaviour.In fact, the presence of a minimum and a maximum radius for the growth of the particles acts as a limiting factor, causing the radii of the circles to have the shape of a Weibull distribution (not shown).

Comparison of the growth mechanism with other models
The model of growth introduced in [6] and used in this paper belongs to the field of stochastic geometry, therefore, it is relevant to compare it to similar mathematical methods mentioned in the scientific literature.
The conditions for growth defined in relation (1) and used in [6] follow the so-called touch-and-stop model of growth [2], which can be defined as a pattern formation process intermediate between the random sequential absorption (RSA) model [2,7] and the Johnson-Mehl-Avrami-Kolmogorov (JMAK) model [2].In fact, in the touch-and-stop model of growth, the particles are characterised by shape persistence (as in the RSA model), but not by size persistence (as in the JMAK model) [2].
In addition, the model introduced in [6] adds a new feature to the touch-and-stop model, as it allows new centres of growth to be seeded once all particles have stopped growing.This is done in order to obtain a lower planar void area, as it was found to be necessary for the generation of realistic representations of granular materials [6].
It is worth mentioning another similar growth mechanism called the lilypond model, which is characterised by hard grains growing radially with a constant speed [7,16].This method, is based on a single finite set of starting points in a subset of R 2 or R 3 , thus it does not involve multiple generations of growth as it is done in [6].In addition, the lilypond model possesses the smaller grainneighbour property [7,16], which states that every grain has to touch at least one other grain that has a smaller or equal radius.In the model described in [6], this feature cannot be reproduced, because the algorithm developed aims at reaching a specific planar air void content while respecting the conditions described in relation (1).This aspect can be seen in Fig. 1 and it happens because the execution of the algorithm stops when a chosen planar void area is reached, even if some of the seeded centres have not yet touched any other neighbouring particle.Another reason for this phenomenon is the fact that boundaries are imposed for the growth of particles, thus, if any of the points defining a particle touches an edge of the domain, its growth is stopped, even if the grain is still isolated.This concept is clearly defined by the second condition in relation (1).
As an example of DEM model, the Void Expansion Method (VEM) can be mentioned [21].The VEM follows a growth principle similar to the ones in [2,6], and in the lilypond model, as a number of particles called structural particles are cyclically grown in three dimensions in combination with void particles.The difference, however, lies in the fact that after each step of growth an equilibration of the assembly of particles is executed by the means of chosen physical principles, e.g. a linear elastic contact law between the particles and damping [21].
The main reason that led to the development of the method introduced in [6] and analysed in this paper is the intrinsic complexity of the existing methods, which arises from the fact that they are based either on complex mathematical formulations or on the implementation of physical laws that rule the interaction of the particles in the computer models.

Combination of 2D packed planes to obtain 3D virtual air pores
Each 2D domain created as explained can be converted to 3D, i.e. circular particles can be converted into spheres, and ellipses into ellipsoids, each with its centre at the same place as the original circle or ellipse.The aim of the procedure is, however, not to produce a layer of randomly distributed spheres or ellipsoids, but to produce a thicker virtual representation of a granular material.Therefore, multiple such layers of spheres (or, in principle, of ellipsoids) can be stacked on one another at fixed spacing in order to obtain a structure resembling that of a granular material.The z coordinate of each x y plane containing the centres of the virtual particles, k, is set in the automation of the stacking process as: where h max i is the radius of the largest sphere in the ith plane and the value of μ has to be found based on the properties of the material that needs to be generated.In the case of asphalt, the value of μ was set as 0.8 and it was found with an optimisation algorithm whose objective was to generate a realistic structure for the void space in terms of voids size and isotropy.The result of the stacking process is that many particles belonging to the different planes are intersected and between them a portion of empty space is left.The stacking process described so far is a pragmatic method to obtain realistic representations of the air voids in granular materials using only geometric principles.Therefore, it is not meant to represent any of the natural phenomena behind the actual manufacture of such materials, e.g.friction or contact forces.The validity of the approach is discussed in the next sections and its use is motivated only by its effectiveness.
In this paper, we consider asphalt as a granular material in order to perform a comparison with real samples.The size of the particles and of the planes, however, can be adapted (reduced or increased) to represent different granular structures.
The use of spheres to describe asphalt may seem to be unrealistic, as the aggregate composing the material almost never has such a regular shape.The forced overlapping of the spheres, however, means that the intersected particles become part of the portion of space representing matter in the virtual asphalt mixture and, consequently, the solid fraction is a highly non-spherical agglomeration with a surface comprised of multiple spherical caps.In this paper, the aggregate is not analysed, as only the void space is studied.However, the intersection of spheres or ellipsoids belonging to different planes is meant to generate a solid structure that resembles a packing of real stones.
It is relevant to add that from the point of view of stochastic geometry it is possible to consider the result of the process just explained as a germ-grain or Boolean model Ξ [18] defined as: where Ξ n are the compact subsets of R 3 defined by the spheres (or ellipsoids) of radii r n obtained following the growth mechanism described in Sect.2.1 and x n are their centres [7].For a deeper discussion of this aspect, see [7,18,19], and [22].
For convenience, the method described so far will be referred to as the Intersected Stacked Air voids method or ISA method.

Locating the void space
In the approach explained in Sect.2, the intersected particles cannot be considered singularly as stones (as done in DEM models [4,9]), because their original shape is lost, but their surrounding void space may be more representative.
For this reason, an algorithm was developed to locate the void space inside a virtual granular material sample.The void location algorithm starts by generating a set of points that covers the whole volume occupied by the virtual sample.These points are compared with the particles associated to matter in the virtual sample and converted into a 3D Boolean matrix that has zeros in the void space and ones in the matter space.The sampling process allows the generation of a 3D map of the void space, as seen in Fig. 2 with a much exaggerated vertical scale.In Fig. 2, the distance between the sampling planes on the z axis is 0.4 mm, i.e. the same distance used to obtain the X-ray CT scans used in this paper.The sampling grid that is used shows numerically what could be deduced by visual inspection of the surface mesh mentioned above, i.e. that the very first and last layers of the virtual material need to be discarded from the analysis, because they contain particles that are "floating" below or above the solid matter due to the mechanism of the algorithm.
In order to show the effectiveness of the sampling method used combined with the algorithm described in [6], a comparison between the CT scan of a real sample and a slice of a virtual sample is shown in Fig. 3.A simple visual inspection of the real and the virtual CT scans shows a high degree of similarity in the shape and size of the void areas in both images (black areas in Fig. 3).In addition, the void patterns in the virtual slice look compatible with the ones seen in the real slice.The visual inspection, however, must be followed by the more strict analysis done in the next sections.

Grid spacing the sampling matrix
The distance between the points in the grid imposed on the 3D virtual model can be changed and a parametric analysis can be performed.The grid spacing, ε, was analysed for the values where D domain is the diameter of the 2D domains used for the generation of the 3D model (D domain = 10 cm in this paper).The lower boundary for ε was found as the limiting value that could be used on a normal office computer with a computational time lower than 45 min.The analysis performed for a virtual sample with air void content equal to about 9 % is shown in Fig. 4. The data represented in Fig. 4 shows that a small difference exists between the planar air void contents evaluated with a variation of the grid spacing.In particular, for the curves shown in Fig. 4 an average air void content of 8.75 % with standard deviation of 0.52 % was found.Therefore, it is possible to state that the result of the sampling algorithm is grid independent, thus, the code can be run according to the specific needs of the user.In fact, if a very detailed grid is needed the user can choose a small value of ε, while if a lower number of points is sufficient ε can be set to a higher value.The use of a dense grid also allows a more precise analysis of the generated 3D domain, as some portions of the void space may be too small to be effectively located with a high spacing between the sampling points.
Generally speaking, in the case of a circular domain the number of sampling points in the x and y directions can be found as The number of sampling planes in the z direction, n z , is also a function of the value of ε and it depends on the thickness of the virtual sample.Finally, it is relevant to notice that the set of points describing the air void content from Fig. 4 is in the shape of a so-called "bathtub" curve [17].This behaviour is a first validation of our approach, since it was found in both real and virtual samples of granular materials used in civil engineering [5,17,25,26].

Topology of the void space
The sampling matrix used to discern the void space from the matter in the virtual models is a cubic lattice.Therefore, each point that is not on a border of the lattice will be connected to 26 other points [27].The number of connections of each point in the porous space to other neighbouring points belonging to the porous space is called the void coordination number, C n [24].As reported in [27], even if the maximum void coordination number is 26, real porous media have a much lower value of C n .For example, in [24], average values of C n between 2.90 and 10.64 are reported and successfully compared to real values from the literature.
The study of the coordination number is of interest in this study, because it determines how neighbouring points belonging to the pore space are connected to one another.The behaviour of the void coordination number for a dense and a porous samples generated with the ISA method and sampled as described in Sect. 3 is shown in Fig. 5.For comparison purposes, the coordination number of real samples is also shown if Fig. 5.The data shown is satisfactory, as in the virtual dense sample (9 % air void content) most particles have a low C n , while in the virtual porous sample (25 % air void content) many particles have higher values of the void coordination number.In addition, the analysis of Fig. 5 clearly shows that the virtual materials have a realistic coordination number when compared to real samples.The data gathered in Fig. 5 also shows the variability of the results Fig. 5 Void coordination number versus normalised number of points with a given coordination number 123 obtained from different simulations.In particular, 15 virtual models with 9 % void content and 15 virtual models with 25 % air void content were generated and the markers in Fig. 5 represent the average values obtained for each coordination number.
In order to show the general validity of the ISA method the number of points in the void space with a given coordination is normalised, thus, allowing the comparison between arbitrarily thick samples.The normalisation is also necessary because the CT scanned samples come with a given resolution that is determined by the chosen equipment, thus, it is usually not possible to match exactly the number of pixels in the real samples in the virtual models that are generated: this mismatch would lead to curves that look different because the number of points under analysis is different.The normalisation is achieved by dividing the number of points with a given coordination number by the total number of points belonging to the void space.Let us specify that to achieve a successful comparison it is necessary to use a small value of ε, which allows a very good resolution in the virtual slices seen e.g. in Fig. 2. In fact, while for the calculation of the air void content the resolution is not a concern, it is when considering more complex properties of the material.As a rule of thumb, values of ε equal to about 0.5 % of the longest side of the domain were found to be effective.
In addition, it is worth mentioning that if the coordination number of the void points is computed for selected slices inside the material on the z axis it will have a very low variance (not shown).
Finally, it should be noted that the void coordination number may show some degree of variability, due to the fact that the virtual and real specimens never show the exact same characteristics and to the techniques needed to threshold the CT scans to isolate the void space.
Additional data concerning the void coordination number of the void space in the virtual samples built for the present study are shown in Table 1.The comparison between the values of the average void coordination number, C n , with those from [24] shows that the void coordination numbers obtained with the ISA method are realistic and similar to the results obtained in other studies.Since the present work mostly concerns the generation of void channels, a successful comparison of the void coordination number is a good sign that the model is able to produce realistic results.This measure, however, does not provide information on the pore network as a whole, as it considers points with maximum distance of √ 3 ε due to the choice of the cubic lattice.For this reason, further characterisation methods are shown in Sect. 5.

Analysis and comparison with real asphalt samples
The results of the model developed are now compared with real asphalt samples.To begin with, it is relevant to notice that the virtual samples generated for the present study show a realistic behaviour in terms of internal structure.In fact, not only the average air void content in the samples follows a "bathtub" shaped curve as seen in Sect.3.1, but also the air void content in the radial direction shows the expected behaviour.In fact, as shown in Fig. 6, the air void content in the centre has a value that is very similar to the average air void content of the whole sample, while it increases significantly approaching the border of the domain, as it was also observed in previous studies of real and virtual samples [5,26].

Isotropy in the virtual samples
In the previous studies using virtual granular materials [11,24], isotropic media were considered.Since our aim was to compare the virtual samples to asphalt samples, we developed an anisotropic 3D model.In fact, asphalt is isotropic on the x y plane, but anisotropic if the x or y directions are com-  pared to the z direction (i.e. it is a cross-anisotropic material), as shown from studies on its hydraulic conductivity [14] or mechanical testing [13].
The number of connections between points of the sampling matrix in the void space, δ, can be calculated for all the virtual domains generated.The number of connections in each direction was divided by the total number of points in the void space in order to obtain a normalised value.This was done as in the 3D Boolean matrix generated in the analysis of the void space the number of sampling points in the z direction is generally different from the number of points in the x and y directions.A normalisation of the number of connections is interesting also because it provides a measure of the global connectivity of the porous space [24].This measure will generally depend on the chosen grid spacing ε, however consistent comparisons can be performed if the value of ε is fixed.
Based on the parameter δ, the isotropy of the virtual domains can be assessed.A simple comparison of the results gathered in Table 2 shows that between δ x and δ y the relative difference is lower than 0.3 %, while the relative differences with δ z are up to almost 2.8 %.Therefore, anisotropy on the z axis only is shown in all the virtual samples.
The analysis of the data generated for this paper strongly suggests that by purposely changing the value of μ (Eq. 3) the isotropy of the virtual materials can be controlled and adapted to the specific needs of the investigation.Future research should investigate the effect of the variation of μ on the properties of the 3D models that are generated.
If the model does not need to show anisotropy in the z direction, a different approach can be used.The three dimensional model could be generated directly rather than by stacking planes of extruded particles by following a procedure similar to the one described in Sect. 2.

Statistical analysis of the size of the air voids
In order to analyse the size of the air voids two approaches can be followed: the first one is the use of the maximal ball (MB) algorithm mentioned in [24] and involves the growth of particles in the void space in order to estimate their maximum size and their coordination number, while the second one Fig.7 Reconstructed air void channels in a virtual sample consists in the actual reconstruction of the void channels.Since our aim was the generation of void patterns in porous media and we already determined the coordination number, we used the second approach.An example of the result of the reconstruction of the void channels can be seen in Fig. 7, where the void space in a virtual sample is shown (a region of interest was selected for clarity).In Fig. 7, void channels connect the top of the sample to the bottom, while in cases with lower values of air void content the reconstruction only consists of a number of unconnected air voids of various sizes in the domain (not shown).
Moving from the sampling matrix described in Sect. 3 to an actual reconstruction of the void space allows a clear comparison with real samples.
We define the void volume, V void , as the volume in voxels of each portion of void space that can be considered as a single entity.In particular, two neighbouring portions of space (subsets of R 3 ) are considered connected only if they share at least a face.
Since the samples have different values of air void content, and, therefore, different void sizes, we compared their voids based on a dimensionless parameter defined as: where V max is the volume of the largest void in each sample.
After evaluating the volume of the voids and the variable λ for the reconstructed models, we were able to establish that the relative void size in all our samples can be described with a Weibull distribution, as shown in Fig. 8 (sample with 9.2 % air void content).
In fact, the approximately linear behaviour of the data shown in Fig. 8 suggests that a Weibull distribution could be used to fit the empirical values.If the same data used in Fig. 8 were to be used in their unnormalised form, they would still show the same approximately linear behaviour in the Weibull plot, since the normalisation was introduced only to visualise dimensionless results.
In [17], the authors analysed the X-Ray CT scans of real asphalt samples and showed that in this material the void size distribution can be described with the a Weibull model.Therefore, recalling the notation used in [17], the where λ is the data to fit, β is a shape parameter, and θ is a scale parameter.The corresponding cumulative distribution function (CDF), shown in Fig. 9 for a virtual and two real samples, can be written as: The maximum likelihood (ML) estimates of β in Table 3 are similar (slightly lower) to those reported in [17], but the value of θ cannot be compared, as the curves fit different kinds of data.However, a similar shape parameter β implies that the curves have overall similar shapes, thus, the comparison is acceptable.
In fact, for our model the actual size of the voids (which is described by θ in the Weibull equation) is not of great importance, since the virtual samples can be very easily scaled and adapted, while the ratio λ void is very significant and actually determines the realism of the computational reconstruction of the void space.In other studies such as [24] different probability density functions were used to fit the air void size distribution, however, that was done because different kinds of materials were analysed (e.g.isotropic rather that anisotropic).The goodness of fit of the Weibull distribution for the computed data was evaluated with the Kolmogorov-Smirnov (KS) test [8] at a 5 % significance level (0.05).The KS test accepts the null hypothesis that our empirical results come from the respective Weibull interpolating curves, obtaining p values between 0.14 and 0.63 (>>0.05).
Because the meshed pore space is described digitally it has the potential to be readily converted to a mesh for finite element or finite difference modelling, e.g. for the computational analysis of fluid flow through the pore space.As an example, Fig. 10 shows a virtual reproduction of the operation called asphalt coring obtained using the ISA method.A visual inspection suggests that the virtual model looks realistic, even if the single stones are not represented.It should be noted that what Fig. 10 shows is not the void volume as seen in the other figures in this paper, but the portion of space associated to matter in the 3D model.This was done to show that further developments of the ISA method will lead to the study of a porous material as whole, thus, also taking into consideration the aggregates that are here neglected.

Analysis the specific surface area of the reconstructed virtual samples
From the point of view of stochastic geometry, it is possible to analyse a subset W ⊂ R 3 of a random model Ξ (see Sect. 2.3) to define its basic characteristics, e.g. the volume density, V V , and the specific surface, S V .The subscript V is used to indicate that the analysis is based on a chosen sampling window corresponding to the subset W and with volume V (W ).As reported in [18], the volume of the portion of Ξ lying in the observation window W is a random variable with expectation EV (Ξ ∩ W ) = V V V (W ), where V (Ξ ∩ W ) is the volume of Ξ restricted to the chosen observation window.In a similar way, the specific surface (or surface density) can be defined as the density of the random surface measure [18].In particular, if S(Ξ ∩ W ) is the surface area of Ξ restricted to W it is possible to write that ES(Ξ ∩ W ) = S V V (W ).
In this paper, the analysis is focused on the porous space, therefore, these properties are evaluated for the air voids.In particular, when the observation window W is the whole sample the volume fraction is dimensionless and represents the air void content of the virtual or real specimen.The calculation of the specific surface area is also based on the air voids and can be performed either by using 3D models directly or by using mathematical methods such as those explained in [18].Since the aim of this study is to provide a first validation of the ISA method, the value of S V is here calculated by elaborating a number of 3D virtual models as this is theoretically simpler.For the analysis of the aspects described above, five 3D models for five random air void contents in the interval from 5 to 30 % were generated in order not to influence the comparison.The values shown in the figures described below represent the average of the results obtained for each value of air void content.
The results of the calculation of the specific surface area are gathered in Fig. 11, where it is shown that the virtual and real samples lie close to one another and have a linear behaviour (the observation window is the cylinder corresponding to the whole sample).The reliability of this calculation, however, is limited to the X-ray CT scans of the real samples available to the authors, thus, further investigations should be performed to check how well this method is able to reproduce different kinds of asphalt mixtures.
In the field of construction materials, a slightly different method is generally used to describe the properties of air voids systems, as the void surface area, S(Ξ ∩ W ), is divided by the total air volume, V (Ξ ), rather than by the volume of the observation window W [1,15].For completeness this parameter was also calculated and shown in Fig. 12.The The data shown in Figs.11 and 12 allows a further discussion of the behaviour of air voids in granular materials.In fact, it can be observed that the surface of the voids increases with the air voids content when compared to the total volume, however, when this parameter is compared to the volume of the voids it shows a decreasing trend.This decreasing trend is not caused by decreasing values of the surface or the volume of the voids, but it is the result of two reasons: (a) large objects (in this case, voids) always have smaller specific surface areas (meant as the ratio between area and volume of voids) due to the reduced number of objects in a certain volume, and (b) larger voids are more interconnected so there are less surfaces in such voids.Generally speaking, it can be observed that high values of specific surface area defined as in Fig. 12 correspond to fine air void systems [1,15].
Finally, it is relevant to add that the specific surface is not able to describe the size distribution of the air voids, i.e. it cannot describe the number of void particles with a given Fig. 13 Distance of points in the matter space from the closest boundary volume it is a single index used by the industry to quickly characterise the properties of multi-size materials [15].However, this parameter can be used as an indicator of the air void distribution for void systems with a similar total air void content [1].

Analysis of the matter space in real and virtual samples
As the observations about Fig. 10 made in Sect.5.2 are only qualitative, it is interesting to calculate some quantitative data to perform a preliminary characterisation of the portion of space associated to matter in the 3D models.For this purpose the Euclidean distance transform [19], EDT, was used to analyse the matter space.For each matter point the EDT yields the distance to the nearest void point.Therefore, the EDT was used to assess if the distance of matter from the voids in the generated models was compatible to the same metric applied to real samples.This approach was applied to all points in every layer of the virtual models, thus, the results do not depend on additional assumptions made by the authors.In Fig. 13, the maximum and average distances of the "matter" points from the void boundaries are represented for a number of real and virtual samples.By observing the data shown in Fig. 13 it can be concluded that the virtual samples fall on the curve interpolating the results for the real samples with very small deviations.The fitting curves calculated from the data from real samples were found by using a power law with two parameters (robust least-squares regression, R 2 = 0.9961 and R 2 = 0.9789 for the average and maximum curves, respectively).
It is relevant to add that in order to allow this kind of comparison the input files or binary images need to have the same resolution, otherwise a relative indicator (e.g.distance over maximum distance) should be used.

Conclusions
In this paper we presented the analysis of a stochastic model based on geometric principles that can be used to build virtual representations of the porous space of granular material.The following conclusions can be drawn: -The Intersected Stacked Air voids (ISA) method can be used to generate representations of the porous space of granular materials with realistic features.-The comparison, qualitative and quantitative, with real samples provides validation of the approach analysed in this paper.-Both continuous and isolated voids can be generated the latter becoming the most common as the void ratio decreases, as found in genuine samples of asphalt.-The void size distribution of the anisotropic virtual samples can be modeled with a Weibull probability distribution.-Meshed versions of the virtual porous space can be created for computational modelling, e.g.computational fluid-dynamics (CFD).

Fig. 2
Fig. 2 Sample layers obtained with the sampling algorithm (porous space is black, 0.4mm vertical spacing)

Fig. 3
Fig. 3 Comparison between a real and a virtual CT scans (porous space is black).a Real sample, 20 % air void content, 100 mm diameter.b Virtual sample, 21 % air void content, 100 mm diameter

Fig. 6
Fig. 6 Sampling planes in a typical cross section (a) and example of radial air void content in different sampling planes (b)

Fig. 11 Fig. 12
Fig. 11 Surface of void space divided by volume of the observation window versus air void content

Table 1
Average void coordination number for the generated virtual samples, C n

Table 2
Normalised number of connections, δ (directions from Fig.

Table 3
ML estimates of the parameters for the Weibull model (Eq.8) Fig. 10 Virtual reproduction of asphalt coring (100 mm diameter)