1 Introduction

In order to design civil engineering materials, such as asphalt, with optimised workability, durability and mechanical properties, it is of the foremost importance to understand how size distribution and shape of aggregates influence their porosity, stone interlock and amounts of asphalt binder or mortar required for an effective binding action [1, 2].

In order achieve a proper design of aggregate-based materials, the specifications allow ranges of gradations and shapes of the particles. There is no standard to determine how particle changes influence the overall material behaviour [3]. This is because laboratory studies of material behaviour are expensive and time consuming. On the other hand, most of the computer simulations of granular materials at the scale levels, still have drawbacks such as generating realistic randomized 3D particles, which are difficult to generate [4,5,6] and the computational cost involved in quantifying the interaction between numerous particles [7, 8].

The first step to create computer simulations of aggregate-based materials is to reproduce realistic aggregates virtually. For example, [9] employed a spherical harmonic function to reproduce over ten thousand of real-shaped aggregates based on the reconstructed aggregate data by X-ray CT technology, and further quantified their shape and size [10, 11]. [13] created realistic particles by assembling spheres, which is could be very costly computationally. [14] created angular particles which, although appear realistic do not have morphological properties, such as angularity or sphericity controlled. [15] simulated aggregates using revolution solids. [16] deformed aggregate projections obtained from a CT Scan gradually, until matching with a target distribution. Garcia et al. similarly, [13, 14, 17] were able to generate aggregates with concavities by creating the particles based on surface coordinates generated by a spherical harmonic series.

After the virtual aggregates are available, these need to be packed. A very extended methodology to pack virtual aggregates is the Discrete Element Modelling (DEM) method, which considers a granular material’s discrete nature, has limitations when it comes to incorporating particles with realistic shapes. The reason is because it is computationally expensive, since a granular assemblage is composed of a multitude of particles interacting [5]. Many studies have modelled packings of aggregates with realistic shapes, using different Discrete Element Methods [19,20,21]. However, the computational cost of DEM is also very high. For that reason, these authors propose to develop new techniques that allow a realistic modelling of particle size and morphology, for a better and deeper study on how they influence different materials behaviour.

An alternative to DEM methods are Physics Engines, which are mainly used in the videogame and animation industry. These serve to perform simulations based on real physics parameters. Examples of previous related works are those from Peng Cao et al., who modelled convex polyhedrons as coarse aggregates [6] and created their convex hull using the Quickhull algorithm [7]. Aggregate conglomerates were created by importing the particles to the Physic Engine SWIG-Bullet multi-body, based on the Monte Carlo method [8]. Izadi et al. simulate pluviation phenomenon using the Bullet physics engine as well [22]. The gravel was represented through randomly shaped polyhedrons, created with the Voronoi’s tessellation [23]. Thus, it can be concluded that a physics engine can be employed as a geotechnical engineering simulation tool. In addition, Garboczi et al. [12] and Zhu et al. [13] developed algorithms that work similarly to those in a physics motor to realise random pacing of realistically shaped aggregates. From these papers, the sections about the overlapping between particles will be relevant for the adoption of Physical Engines to pack virtual aggregates. Finally, [16, 24] have also used statistically-based methods to generate packs of particles from planar 2D surfaces.

The objective of this paper is to show a novel algorithm for the creation of virtual realistic 3D aggregates. The algorithm uses inputs from 2D aggregates to define the 3D geometries. Previous research has shown that there are linear correlations between 3 and 2D size descriptors [28]. The algorithm consists of a spherical density function and a Perlin noise to deform the spherical density function. By combining both density functions, a scalar field is created, which generates a cloud of points that re-creates the shape of the aggregates in a realistic way. Then, the points are triangulated by using the Marching Cube algorithm [25,26,27]. Finally, the particles are deformed until they match the aspect ratio and minor Feret of the aggregates that we want to deform. The model proposed also offers the possibility of changing the particle sizes and shapes by considering the statistical parameters of real aggregate samples.

2 Materials and methods

2.1 Experimental measurement of aggregates’ morphological properties

Nine types of coarse aggregates from different sources and with different shape characteristics have been used in this study, namely glass spheres, crushed glass, round gravel, limestone of maximum size 14, 10, and 6 mm and two types of granite of 14 mm. Each type of stone had its morphological information extracted in order to implement the probabilistic parameters in the virtual aggregates model. The shape factors used to measure the aggregates are described in Table 1.

Table 1 Morphological parameters to characterize particles
Full size table

Pictures have been obtained from the aggregates resting over their most stable position and, the analysis of their geometries has been carried out by the digital image processing software ImageJ, see examples of the aggregates in Fig. 1. The image analysis included (i) establishing the scale, (ii) highlighting the borders of the particles and correcting the shading following the procedure mentioned in [32] and, (iii) measuring and calculation the properties from Table 1, using the BioVoxel plugin and its functions for extended particle analysis [33]. The distribution of the properties from Table 1 has been analysed in the Minitab software using the Weibull distribution, which is a commonly accepted one to create aggregate distributions. The morphological information of the aggregate types has been quantified in terms of the 50th percentile, P50.

Fig. 1
figure 1

Examples of the aggregates analysed

Full size image

Finally, the height of the aggregates has also to be adjusted to that of the aggregates. It has been defined by a new parameter, \(\beta\), which we will call thickness from now on, and is the ratio between the minor Feret and the average height of 50 stones, measured with a calliper. The height has been measured when the stone is lying in the most stable position, in the direction perpendicular to the surface.

2.2 Physics engine: Unity 3D

The platform chosen to do this research is Unity3D, which utilises NVIDIA´s PhysX physics engine, and C# as the development language. PhysX is a multi-threaded physics simulation, which supports real-time simulation of rigid and soft bodies, fluids and springs, with a focus on application to the gaming industry [34]. This physics engine is well documented and allows to detect collisions between objects with complex shape thanks to its built-in collision detection and response algorithms and mesh collide [35].

2.3 Create complex virtual particles with dents similar to that of aggregates

In order to re-create dented virtual aggregates, this section presents the spherical primitive + Perlin noise method, will deforms a spherical primitive using a Perlin noise matrix, until its angularity and circularity, their denting level, becomes like that of an aggregate. The Perlin noise is widely used in the animation and movies industry in order to simulate natural-looking landscapes. In general terms, it consists of a generator of pseudo-random numbers, which are softened by a weighted average of the points [36]. Figure 2 summarises the process to produce a virtual aggregate.

Fig. 2
figure 2

Process to create a dented particle

Full size image

Firstly, a spherical primitive with radius r, is inscribed in a grid of N x N x N voxels, see Fig. 2a, where N is the number of voxels and d, the distance between the nodes of the voxels. The radius of the sphere is defined by the following equation:

$$r = d \times \frac{N}{2}$$
(4)

An increase in the number of voxels, N, implies an improvement of the resolution of the 3D model but has the inconvenience of a higher computational cost. In this paper, a grid of 7 × 7 × 7 voxels with 1 unit distance between them has been chosen by trial and error, since it has been found that this permits a good resolution while avoiding large memory requirements [37].

Then, a 3D density function \(f(x, y, z)\) that determines which lattice nodes, with coordinates, \(x, y, z\), are inside or outside the sphere is defined. In this case, the nodes outside the sphere are positive, the nodes inside the sphere are negative and, the nodes at the surface have a value of 0. The formula that defines the weights in the density function is (5). In addition, a weight parameter (\({S}_{c}\)), has been added that increases the density values of the spherical functions. This factor is used to control the influence of the spherical density function in the final shape of the aggregates. In other words, it controls the sphericity of the virtual aggregates. As the numbers of voxels is 7, and the distance between them 1, the radius, from the centre of the voxel, is 3 and the density function values inside the sphere range between 0 and 3.

$$f\left( {x,{ }y,{ }z} \right) = \left( { - r + v} \right){*}S_{c}$$
(5)

where \(v\) is the distance from the voxel’s centre to the centre of the sphere.

In order to introduce dents in the particles that are like those in crushed aggregates, a second grid with 7 × 7 × 7 voxels is generated and combined with the spherical one, see Fig. 2b. The node values, with coordinates \(x, y, z\), (\({R\left(x,y,z\right)=P}_{v}\)), of the second grid are obtained from a Perlin noise generator, which values may range between 0 and a maximum value of 4, depending on the degree of deformation desired, we have called the highest value of the Perlin’s noise range, f. The advantage of using a Perlin noise is that it produces a gradient, natural looking surfaces [39]. If the value is 0, there is no deformation. If the value is 4, the deformation is the highest. This value has been determined by the trial and error method and, a higher value could produce unrealistic shapes (see Fig. 3), with holes. Each aggregate has a different matrix. To control the angularity of the aggregates, \({P}_{v}\) can be multiplied by a factor that we have named Intensity of deformation, \({I}_{def}\), which weights \({P}_{v}\).

Fig. 3
figure 3

Example of un-realistic and realistic virtual aggregates

Full size image

Finally, the density function of the aggregates, \(F\left(x,z,y\right)\), is created by multiplying the spherical density function and the noise density function as follows, see also Fig. 2c:

$$F\left( {x,z,y} \right) = f\left( {x,{ }y,{ }z} \right) + R\left( {x,y,z} \right) = \left[ {\left( { - r + v} \right){*}S_{c} } \right]*\left[ {P_{v} {*}I_{def} } \right]$$
(6)

Finally, a Marching Cube algorithm is used in order to triangulate the shape of the virtual saggregate. The points in the voxels’ edges that cross the density values equal to 0 are identified by the Marching Cube and triangles are generated that transverse them and form the shape of the virtual aggregate.

In occasions, the method presented above may produce virtual aggregates with topological inconsistencies, such as holes in the mesh, see Fig. 3 as an example. In this research, it could happen for two reasons: (i) due to the intrinsic nature of the Marching Cube algorithm, which in occasions may interpret that some faces correspond to holes [17] and (ii) the random density matrix which can deform the particles excessively. These ambiguities occur when there is more than one possible way to triangulate the points at the voxels’ corners and these can be interpreted as holes. A solution for this could simply be to increase the number of voxels in the matrix, N; other researchers have also identified this problem and proposed solutions that minimise but, they do not correct it completely [38].

Such virtual stones are discarded by a particle control routine that uses an operation from Unity that “sends a virtual laser ray” that can detect surfaces. In this way the surface of aggregates is mapped and aggregates with holes are detected, depending on whether some rays go through the stone, and when any are found, the stone is discarded.

2.4 Deform the particles to adjust their minor Feret and aspect ratio in the most stable position to that of aggregate projections

The minor Feret and aspect ratio are enough to characterise the dimensions of each particle projection [18]. The algorithm described in the previous section produces angular particles which minor Ferets and aspect ratios do not fit those of the actual aggregates. For that reason, a method has been created that rescales the angular particles, to achieve the desired minor Feret and aspect ratio.

As it will be shown in the Results section, the aspect ratio and minor Feret of each aggregate type have been characterised using the Weibull distribution function, as it was done in [18]. This distribution is defined by two parameters, one of shape, \(k\), and another one of scale, \(\lambda\), which have been measured experimentally, see also Table 2.

Table 2 Input parameters to control size and shape of virtual aggregates in the model
Full size table

To select the aspect ratio and minor Feret of each particle, a pseudorandom function, with range between 0.05 and 0.95 generates numbers, which are imputed as probabilities in Weibull distribution functions of aspect ratio and minor Feret that uses as scale and shape factors measured experimentally. The results are values of minor Feret an aspect ratio which fall within the same distribution of aggregates measured experimentally.

A virtual laser is used to scan every virtual particle and select the different views of the particle, see Fig. 4, based on the position of the points on the particle’s surface “hit” by the ray. Then, an array with the coordinates of the outline plan view is created, which enables definition of the perimeter. From this array, the minor Feret and the major Feret are identified and measured, see Fig. 4d. The area is measured based on the triangular division formed by the segments between the perimeter array and the centre of the particle’s projection, Fig. 4d, following the procedure described in [40]. Note that these particles are still approximately spherical.

Fig. 4
figure 4

Orthogonal views of virtual aggregates and shape factors

Full size image

Finally, one of the views is selected and, the virtual aggregates are rescaled so that their planar projections have the minor Feret and aspect ratios within the Weibull probability range measured experimentally. In addition, the thickness of the aggregates is adjusted so that their thickness fits that determined experimentally.

2.5 Select input parameters that define the size and shape of virtual aggregates

A summary of all the input parameters required to reproduce the aggregates can be found in Table 2. In order to calculate automatically the input combinations, we have used a DE algorithm [41]. The outputs that the algorithm has used are \({P}_{50}\), which is the median value of the perimeter and \({A}_{50}\), the median value of the area.

2.6 Explanation of the use of virtual aggregates through packing.

G2 materials were used to prepare an aggregate packing to explain the use of the virtual aggregates. The packing weighted approximately 500 ± 5 g. The aggregates were placed in a cylindrical container to compact the aggregates, whose dimensions are 9.4 cm in diameter and 8 cm in height. The compaction test was conducted using a vibratory table with flexible springs and an engine speed of 3600 rpm, making the table vibrated in a horizontal sine wave. The frequency was 60 Hz, and the three axes’ displacement ranged between 0 and 1.58 mm. The stones were firstly vibrated for 2 min and secondly, a metal cylinder, with weight 2.7 kg and 8 cm of diameter was dropped from a 5 cm height above the aggregates and, the vibration continued for five additional minutes. Finally, the pack of G2 aggregates was CT scanned as described in [42], and the images were treated with the image processing software ImageJ following the procedure described in [43].

An equivalent mass of G2 virtual aggregates was compacted in the Physics Engine by pouring them in a cylindrical container with the same dimensions as the real one. The total number of particles in the simulation was defined by the total mass of the sample and the density, as Unity 3D measures each particle's volume. The vibrations over the virtual cylinder mimicked those of the experiment and, compaction was stopped when the height of the virtual specimen reached that of the experiment. Parameters used in the simulation were concave collider, default contact onset 0.0001 m, bounciness 0.315, static friction 0, 30 iterations per frame and a fixed timestep of 0.002.

3 Results and discussion

3.1 Experimental values of the aggregates studied

The minor Feret and Aspect Ratio data obtained from ImageJ have been analysed using the Anderson–Darling adjustment, which indicates the goodness-of-fit to a distribution curve; the closer the result is to zero, the better the fit of the curve to the data points [21]. Based on [18], a two-parameter Weibull distribution functions has been used to fit the minor Feret and Aspect Ratio. An example can be seen in Fig. 5, which shows the cumulative probabilities for the minor Feret of G1 aggregates. The Anderson–Darling test statistic value, A2, is 0.479.

Fig. 5
figure 5

Weibull plot of Minor Feret

Full size image

Table 3 shows the values that characterise the aggregates studied. See for example the glass spheres, which show aspect ratio approximately equal to 1, minor Feret approximately 15 mm, and area and perimeter like those of the circle with equivalent diameter.

Table 3 Experimentally measured shape factors that characterise the aggregates studied
Full size table

The aggregates with the largest shape parameter values are RG and GS, whereas the aggregates with the smallest shape parameter values are S and L3. The aggregates with a more uniform size distribution are L1 and G1, with shape values greater than 10; followed by S, CG and G2, with shape values between 5 and 10. In addition, the aggregates with highest AR are GS, RG, and L2, while the aggregates with lowest AR are CG, L3 and S. Those are also the aggregates with higher and lower angularity and circularity. Furthermore, the thicker aggregates, with β greater than 0.5 are RG, G2 and L1, whilst the thinner aggregates are L3, S, and CG.

Finally, Fig. 6 shows the aggregate classification based on their circularity and sphericity, using the Krumbein chart used for visual evaluation of aggregates [44]. From results from Table 3, it can be concluded that \({R}_{50}\) and \({C}_{50}\) are linearly related:

$$R_{50} = 1.10\beta - 0.12; R^{2} = 0.91$$
(7)
Fig. 6
figure 6

Aggregate chart classification of the aggregates studied

Full size image

In addition, a linear trend can be also observed between these two parameters, \({AR}_{50}\) and \(\beta\), which implies that, for the aggregates used in construction materials such as asphalt or concrete, its value does not need to be measured:

$$AR_{50} = 0.56\beta + 0.42; R^{2} = 0.86$$
(8)

\(\beta\) does not show any clear relationship with \({Min Feret}_{50}\) and \({Max Feret}_{50}\).

3.2 Input parameters to generate the aggregates studied

The input parameters that control the scalar density functions are \({S}_{c}\), \({I}_{def}\), f. Varying these parameters controls the degree of particle deformation and contributes to achieve various levels of angularity and sphericity. Their combination has been found by trial and error through a visual inspection using Krumbein’s chart classification, which can be found in [44]. Figure 7 shows the input values used to find a wide range of morphologies that corresponds to those commonly found in the aggregates used in civil engineering.

Fig. 7
figure 7

Aggregate chart classification of the aggregates studied

Full size image

Figure 8, shows the virtual aggregates in a Krumbein chart, based on the circularity and sphericity of the aggregates. This Figure has been prepared to give the reader an idea of the general shape of the virtual aggregates and the parameters used to create them. The objective is to show that the methodology shown in this paper can be used to re-create any aggregate from those commonly used in civil engineering. It is remarkable that this methodology has allowed to precisely control the angularity and circularity of the virtual aggregates.

Fig. 8
figure 8

Aggregate chart classification of the virtual aggregates studied

Full size image

To compare the real vs virtual aggregates, these have been simply represented one versus the other. Figure 9 shows the distribution of Minor Feret of G2 aggregates as an example. In general, it was found that the main difference corresponds to the extremes of the distribution. This can induce errors in the representations of aggregates, such as aggregates that are too elongated, or bigger than normal. If these aggregates would be used to make virtual asphalt mixtures or other types of aggregate-based materials, the most extreme virtual aggregates should be eliminated in order to produce correct mixtures. Obviously, the simulation shows a high level of similarity because the virtual aggregates have been made to fit the experimental data.

Fig. 9
figure 9

Distribution of Minor Ferets for G2 aggregates, calculated from simulations and experiments

Full size image

The algorithm has been validated by comparing the real and virtual morphological information of each type of aggregate. The validation has been done by a correlation parameter named reliability index, which compare the median of the virtual stone´s samples to the experimental results, in order to know the quality of fit of the virtual aggregate model with the real data; see Eq. (9). A reliability index value equal to one indicates a perfect adjustment.

$${\text{Reliability index = }}\frac{{{\text{P}}_{{\text{50 virtual}}} }}{{{\text{P}}_{{\text{50 real}}} }}$$
(9)

The reliability indexes are shown in Table 4. As the Weibull's distribution parameters are inputs imposed by the user, those generated by the algorithms are identical or very close to them.

Table 4 Reliability indexes resultant from comparing the simulations with experiments
Full size table

It is worth mentioning that the perimeter, area, circularity and roundness show the lowest accuracy. This is due to (i) errors induced from the selection of the high number of variables involved in the creation of particles, which were chosen by the trial and error method; for example, f, \({I}_{def}\), and \({S}_{c}\) and, (ii) the measurements of area and perimeter by ImageJ, which may be influenced by shadowing and thresholding. In the future, an optimisation algorithm can be used to identify the aggregates more precisely.

3.3 Comparison of a simulated aggregate packing to an experiment

In order to give an idea of the capabilities of using Physics Engines to produce aggregate packings, Fig. 10 shows a comparison between a simulation and a packing of aggregates G2, obtained through X-ray CT scans. The goodness of fit of the vertical distribution of pores has been found via the Kolmogorov–Smirnov Test [45], using the critical value for a confidence level at 95% (\({D}_{crit, 0.05}\)), using:

$$D_{crit, 0.05} { } = 1.36 \sqrt {\frac{1}{{n_{x} }} + \frac{1}{{n_{y} }}} ,$$
(10)
Fig. 10
figure 10

Vertical distribution of air voids in a simulated and experimental packs of G2 aggregates

Full size image

where \({n}_{x}\) and \({n}_{y}\) are the population size of the samples.

\({D}_{crit, 0.05}\) = 0.028; if \({D}_{crit, 0.05}\) is smaller than the maximum absolute distance between the expected and observed distribution functions (\({D}_{s}\)), we can conclude that the values from the experiments and simulations are equivalent. In this case \({D}_{s}= 0.013\), hence, we can confirm that the virtual aggregate’s packing has a very close AV distribution function than the experimental one. Finally note that for the same equivalent mass, the experiment comprised 104 particles, while the simulation 110 particles.

4 Conclusions

This article has shown a novel methodology to re-create commonly used aggregates in civil engineering. It is based on a spherical primitive that is deformed by means of a Perlin noise. Then, the shape and size of the virtual aggregates is adjusted by deforming the shapes until they fall in the same statistical range of the real aggregates. The algorithm has been implemented in the NVIDIA PhysX physical engine, with the objective of being able to use these aggregates in future research that will involve creating virtual aggregate gradations.

The following conclusions have been obtained:

  • This methodology allows to re-create the morphologies of aggregates virtually. Based on the experience of the researchers, hundreds of aggregates were generated in a few seconds.

  • This algorithm has allowed to re-create realistic virtual aggregate shapes that include concavities. The authors are convinced that the precision of these virtual aggregates will be enough to analyse the influence of the aggregate shapes on the packing of aggregates.

  • The angularity and sphericity of the particles are controlled by three inputs that weight the spherical primitive and Perlin noise. In the current research, these inputs have been adjusted by trial and error, until the shape of the virtual shapes was found. In future research, this will be updated by using an optimisation algorithm, which will allow the creation of aggregates in a fast and efficient way.

  • The comparison between the virtual and the real geometrical properties’ data indicates that the algorithm to create virtual aggregates can achieve minor and major Ferets like those measured in real aggregate samples. In addition, this methodology allows to create virtual particles with area, perimeter, sphericity, and angularity similar to those of the real materials.

  • In occasions, the methodology proposed in this article produced aggregates that presented unnatural shapes, including holes and bottlenecks. The reasons for this are still unsolved. The unnatural particles appeared very seldomly. However, if hundreds of particles are produced to create an aggregate packing and some unnatural particles appear results could be affected. For this reason, the authors have developed a new methodology to produce virtual aggregates that is presented in [18].

  • It has been found that the methodology allows to re-create full gradations by using the Weibull distribution to decide the minor Feret of the virtual aggregates. However, in future research the extreme aggregates need to be eliminated, to allow re-creating unrealistic aggregates.