Anm: a geometrical model for the composite structure of mortar and concrete using real-shape particles

The composite geometrical structure of mortar composites can be represented by a model consisting of sand embedded in a cement paste matrix and the structure of concrete by gravel embedded in a mortar matrix. Traditionally, spheres have often been used to represent aggregates (sand and gravel), although the accuracy of properties computed for structures using spherical aggregates as inclusions can be limited when the property contrast between aggregate and matrix is large. In this paper, a new geometrical model is described, which can simulate the composite structures of mortar and concrete with real-shape aggregates. The aggregate shapes are either directly or statistically taken from real particles, using a spherical harmonic expansion, where a set of spherical harmonic coefficients, anm, is used to describe the irregular shape. The model name of Anm is taken from this choice of notation. The take-and-place parking method is employed to put multiple irregular particles together within a pre-determined empty container, which becomes a representative volume element. This representative volume element can then be used as input into some kind of computational material model, which uses other numerical techniques such as finite elements to compute properties of the Anm composite structure.


Introduction
Concrete is a composite heterogeneous construction material, which is composed primarily of coarse aggregates (e.g. gravel), fine aggregates (e.g. sand), cement, and water. The cement paste matrix binds the aggregates together and forms a system that is able to carry loads. Concrete without coarse aggregates is called mortar. In a simple composite model for mortar, the particles are sand and the matrix is cement paste, while for concrete, the particles are coarse aggregates, and the mortar serves as the matrix. However, the particle shape characteristics are usually different for sand and coarse aggregates, although in general particle shape, for a given material class, depends only weakly on particle size for blasted and crushed rocks [1][2][3]. But particle shape can be quite different for different mineral classes of fine or coarse aggregates or material that has been prepared in different Z. Qian Materials Innovation Institute M2i and Delft University of Technology, Delft, The Netherlands ways (e.g. crushing vs. naturally rounded). This requires any realistic 3D geometrical model of mortar or concrete to be able to recognize and use various particle shape classes.
There are many existing geometrical models of particles embedded in a matrix [4][5][6][7], most of which employ regular shape particles like spheres, ellipsoids, or multi-faceted polyhedrons. Recently, a method was described for generating packings of realistic particles based on random field theory and Fourier shape descriptors [8][9][10]. This model was oriented more toward particle packings, where particles touch and mechanics is dominated by contact forces. In mortar and concrete, the sand or gravel particles can be thought of as being suspended in a matrix, forming a composite whose mechanics is controlled by stress transfer between matrix and inclusion, not inclusion-inclusion contact forces. The real particle shapes found in construction of naturally-rounded or blasted and crushed sands and gravels sometimes play an essential role in the accuracy of the model, for computing fracture paths and when the property contrast between aggregate and matrix is high [11]. It has been shown that spherical harmonic series are an effective mathematical tool to characterize the shape of particles analytically, and the procedures to retrieve particle shape characterizations for a given class of aggregates from X-ray CT (X-ray computed tomography) scans have been established [12]. Statistical methods have also been developed to generate new particles based on statistics that have been obtained from a real particle dataset or elsewhere [13,14].
What else is needed to build a multi-aggregate 3D model of mortar or concrete with real particles? One must be able to randomly place multiple irregular shape particles into a predefined container, without overlapping. This procedure has been denoted random sequential packing, random dense parking, or random jamming [4,15,16], and we will use parking or placing interchangeably below. In this paper, a geometrical model with irregularly shaped particles is proposed, which has these characteristics, and is denoted the Anm model [17]. This name is taken from the choice of notation for the spherical harmonic expansion coefficients, a nm , which are used to describe particle shape. The Anm model has been written in the computer language C ?? . Simple applications of the Anm model are made to demonstrate its capabilities for simulating the geometrical structure of mortar and concrete, using both periodic and non-periodic boundary conditions. The multiple inclusion models can also serve as a good starting point for discrete element modeling (DEM) [18] and for finite element fracture models [19].
2 Description and properties of an individual particle with irregular shape 2.1 Mathematical representation of a 3D irregular shape An arbitrary 2D single-valued closed surface of a 3D particle can be represented by a function r(h,/) in a 3D spherical polar coordinate system. The origin of the local coordinate system is placed at an interior point, in the case of particles usually the particle center of mass. Single-valued means that a ray from the origin to infinity will only intersect the surface once, which is equivalent to the definition of a star-like particle [20]. The function r(h,/) is, in principle, impossible to express explicitly for random shapes, but it can always be approximated by a summation of spherical harmonic functions for star-like shapes [21][22][23]. It is clear that almost of all sand and gravel used for concrete have star-like shapes (i.e., are star-shaped) [12,24,25]. The spherical harmonic expansion for r(h,/) is given in Eq. (1), where 0 B / B 2p and 0 B h B p).
Y nm (h,/) is the spherical harmonic function with indices n and m and is given by Eq. (2), with -n B m B n, P mn (cos h) is the associated Legendre polynomial, and i is the square root of -1 [26]: Other versions of the spherical harmonic functions are used in various disciplines, with slightly different conventions, as well as real versions, but this is the definition and convention used herein.
The particle shape is uniquely defined by the set of spherical harmonic coefficients a nm , which in general are complex numbers. The total number of the a nm coefficients up to the degree n is (n ? 1) 2 . There should be sufficient coefficients involved in the spherical harmonic expansion to describe a shape accurately. Of course Eq. (1) is only exact when n becomes infinite. The number of the coefficients that can accurately be used is dependent on the resolution of the original experimental shape characterization. As a general guideline, the degree n should be at least 14 and is usually no larger than 26 for common aggregates in mortar and concrete. The number of coefficients needed for accuracy and how this accuracy is determined is discussed more fully in [12,24]. In this paper, we assume that we have spherical harmonic representations of aggregates that are sufficiently accurate for our purposes. An example of an irregularly-shaped particle, which has been described in terms of spherical harmonic coefficients, is shown in Fig. 1. Particle shape databases have been created for varying classes of powders and aggregates withthe procedures proposed in [12,24,25], either from direct X-ray computed tomographic measurements of particles or from a statistical procedure for generating new particles based on such a measured collection [13,14]. A selection from these databases, contained in the Anm model, so far includes two crushed coarse aggregates (gravel), two crushed and naturally rounded fine aggregates (sand), and one Portland cement. Cements are included since cement particles suspended in water can form an accurate geometrical model of unhydrated cement paste, although the focus of the model and of this paper has been on mortar and concrete applications.

Particle geometric properties and operations
Many geometric properties of the particle shape can be computed once the spherical harmonic coefficients are known, including the particle volume, surface area, integrated mean curvature, moment of inertia tensor, length, width, and thickness [25,27]. Length is the longest surface-surface distance in the particle, width is the longest surface-surface distance in the particle such that width is perpendicular to length, and thickness is the longest surface-surface distance in the particle such that thickness is perpendicular to both length and width [25]. A particle's size can be changed but its shape preserved by multiplying each a nm coefficient with the same scaling factor. The scaling factor s is defined by s = W a /W b , where W a is the particle width after scaling, and W b is the width before scaling. The particle width is preferentially used for computational sieve analysis [12] since it is thought to match best with the usual standard experimental sieve classification of particles [28,29]. The particle volume will of course scale by a factor of s 3 , and the surface area by a factor of s 2 . In the Anm model, rescaling of particle size is a core operation of the overall algorithm.

Parking procedures and algorithms
In the simplest sense, the Anm geometrical model simply places or parks [4,15,16] particles in a matrix. The matrix is considered to be a continuum: mortar surrounding coarse aggregates or cement paste surrounding fine aggregates. An empty container is created to represent a specimen, and then particles are placed one after another into this container, from the larger ones to the smaller ones, following a particle size distribution that is given in terms of a sieve analysis, selecting particles from a particle shape database. A preliminary database, containing several sand and gravel types, has been created from previous data as part of the Anm model. This is not the equivalent of equilibrium packing, but when there is a wide particle size distribution, equilibrium packing and non-equilibrium placing produce similar arrangements of particles [30]. The sieve analysis is usually taken from experiment or some kind of standard for concrete materials. A pseudo-random number generator (the default function included in C ?? ) is initialized with a seed value that is different for each run of the Anm model, thus guaranteeing different results for Fig. 1 An irregular shape sand particle described by spherical harmonics each run. The same pseudo-random number generator was used for all random numbers required for the functioning of the Anm model. Using a different random number generator will be considered for future versions of Anm, since such a large number of random numbers can be used for each Anm run.
The Anm model starts with the largest particles as it would be more difficult to place them if they were processed at a later stage. Each new particle that is placed must not overlap any other particles that have already been placed, since that would be unphysical. Particles are included in a sieve range according to their size, which is taken to be the particle width (previously computed). The widths within a sieve range are assumed to be distributed evenly, with equal probability for any width in the sieve range. A width within the sieve range being processed is picked randomly, using the same pseudo-random number generator that gives a random number between 0 and 1, with equal probability. The random number is scaled to the sieve range limits, determining a random width. A particle is chosen randomly from the appropriate particle shape database, which consists of a collection of sets of spherical harmonic coefficients representing a collection of particles. All particles in the database have equal probability and one is chosen randomly using the pseudo-random number generator in a way similar to how the width was chosen. If there are M particles, then the output of the pseudo-random number generator is scaled to pick an integer between 1 and M in order to choose a specific particle. Note that for the particle size ranges considered in the Anm model, particle shape does not appear to depend on particle size, so that any particle in the shape database can be randomly picked and scaled to fall into a given sieve size range [31].
When placing a particle, to simulate mixing we also need to randomly rotate each particle before attempting to place it at a random location. The operation of a rotation on the particle in the local coordinate system can be done by transforming the spherical harmonic expansion coefficients [22] according to an arbitrary orientation defined by the three Euler angles a, b, and c [26]. Each Euler angle is assumed to be uniformly distributed in its range, and a selection is made of the three angles using three random numbers scaled to the appropriate angular range.
All points within the unit cell are equally probable, so the pseudo-random number generator is used again to choose values of x, y, and z that fall within the bounds of the unit cell. The chosen point is used for the origin of the particle being placed. This origin is always the center of mass of the particle [12]. The chosen particle is checked for overlap against the previously placed particles within a certain distance. This distance is controlled by a local bin structure. The unit cell is divided up into bins, typically 10 3 = 1,000, and the information about which bins are touched by a given particle is recorded. Only those particles lying in bins that have a possibility of overlap with the trial particle are checked. This saves considerable run time, especially once a significant number of particles have been placed. If periodic boundary conditions are being used, then the periodic ghost particles (see Sect. 3.3) are also checked for overlap with existing particles. If no overlap is detected, then the particle enters the simulation box successfully, otherwise the program will try to place it at a new randomly chosen location. The reassignment of the location is subject to a predefined maximum number of attempts. After the number of consecutive failures reaches the maximum limit, the particle will be resized to another randomly selected width within the current sieve range, and the program will try again to place the particle at a new random location. Rescaling of the particle size is also subject to a pre-defined maximum number of attempts. If the size scaling does not result in particle placement, then the particle will be rotated randomly to have a new orientation. If the particle still cannot be placed, then a new shape will be chosen from the particle shape database and the process restarts. If enough particles in that sieve range have been placed, satisfying the pre-determined volume fraction for that sieve, then the next smallest sieve class is processed. However, if after a number of such unsuccessful tries, again a predefined parameter, the next sieve range will be processed even if the volume placed is less than what the user desired for this sieve range.

Overlap algorithm
The key algorithm required in the particle placing procedures is to determine whether two particles overlap. This can be examined by formulating and solving contact equations. As shown in Fig. 2, the mass centers of two particles are O 1 (x 1 ; y 1 ; z 1 ) and O 2 (x 2 ; y 2 ; z 2 ), in a global coordinate system O(x; y; z). Two local coordinate systems are also defined and their origins are placed at O 1 and O 2 , respectively. It is assumed that there is a contact point C(x c ; y c ; z c ), and its local coordinates are (h c1 ; / c1 ) and (h c2 ; / c2 ) in the corresponding local coordinate system. In the particle database, the coordinates, relative to the center of a particle, of the corners of a box that just contains the particle are stored. These are updated when the particle is rescaled or rotated. If, for two particles, these boxes do not overlap, then there is no way that the particles themselves can overlap and no further checks are done. In this way, many of the potential particle overlaps can be handled simply.
The contact point C is located on the surface of particle O 1 , with coordinates defined using the spherical polar angles and origin associated with particle 1, so the following Eq. (3) should be satisfied, if the contact point C exists: The contact point C is also located on the surface of particle O 2 , with coordinates defined using the spherical polar angles and origin associated with particle 2, so Eq. (4) must also be valid, Equations (3) and (4) can be equated and rearranged, which yields Eqs. (5): x 2 À x 1 þ r 2 ðh c2 ; / c2 Þ sin h c2 cos / c2 À r 1 ðh c1 ; / c1 Þ sin h c1 cos / c1 ¼ 0 Equations (5) have four unknowns h c1 ; / c1 ; h c2 ; / c2 , which makes them unsolvable in principle. If these equations were expressed in terms of the three unknown Cartesian coordinates of the contact point, then there would of course be only three unknowns and three equations and therefore the contact equations would potentially be solvable if the contact point C existed. However, it is much easier to express the various derivatives needed for application of the Newton-Raphson method in terms of spherical polar coordinates, hence the four unknowns. The strategy to overcome this difficulty is to fix one of the four unknowns (e.g., / c2 ) at pre-selected values in the range of [0; 2p), and then try to directly solve the three equations with the three unknowns h c1 ; / c1 ; h c2 using the Newton-Raphson method. If a solution is found, then it demonstrates that the contact point exists, and the two particles overlap. If no solution can be obtained for every selected value of / c2 , then it is assumed that the contact point does not exist, and the two particles do not overlap. It is recommended to take at least 20 different values of / c2 evenly distributed in the range [0; 2p). The number of such points is another user-selected parameter of the algorithm. This contact function has been tested against direct observation using VRML (virtual reality modeling language) images and has been found to be accurate for particle surfaces that approach to within a small fraction of either particle radius, approximately 1 %. For particle surfaces that approach closer than that distance, it is unknown how well the algorithm discriminates between overlap and non-overlap.
A more accurate and faster contact algorithm for star-shaped particles, which has been more rigorously tested, is described in Ref. [32]. This contact function will be introduced into the Anm model in its next version [33]. Another potential contact function, based on that described in Ref. [32] exists, but since it uses a numerical triangulation of particle surfaces, it can be applied also to non-star-shaped particles [34]. Since the Anm model for the foreseeable future will be Fig. 2 2D schematic of two random-shape particles overlapping restricted to star-shaped particles, the next version [33] will use the contact function of Ref. [32].

Periodic and non-periodic material boundaries
In mathematical models and computer simulations, the behavior of a small-volume model can often be dominated by surface effects, since the surface area to volume ratio increases as the volume decreases. Periodic boundary conditions are a set of boundary conditions that are often used to make a small model seem larger by removing the dominance of the surface since periodic boundary conditions eliminate the surface. In the Anm model, both periodic and nonperiodic boundary conditions can be used. The periodic boundary permits a particle to pass through the surface of the simulation box and the part outside the simulation box is put on the opposite surface by placing a duplicate particle (ghost particle) with the same orientation, in the periodic position, while nonperiodic boundary conditions do not allow a particle to pass through the surface of the simulation box, as shown in Fig. 3. In the specimen with periodic boundary conditions, if a particle passes through a surface but not an edge or corner, then it generates one ghost particle; if a particle goes through an edge but not a corner, then three ghost particles are created; if a particle passes through at a corner, then the number of its ghost particles is seven. In that way, every part of every particle appears inside the unit cell only once. A combination of periodic and non-periodic boundary conditions, applied to different faces of the unit cell, can also be employed.

Demonstrations of Anm for mortar and concrete
For now, the Anm model is designed for two applications: simulating the arrangement of coarse aggregates (gravel) in a mortar matrix, and simulating the arrangement of fine aggregates (sand) in a cement paste matrix. In either case, the volume fraction of inclusions is limited to less than about 50 %. It should be possible, using the Anm model, to simultaneously place sand and gravel to achieve a total inclusion volume fraction of 60-80 %. However, for this wide size range of particles, approximately 0.5-12 mm, a very large number of particles, probably well over 100,000, would be necessary. The run time needed for the current version of Anm to simulate this situation would be prohibitive. A newer, faster version of Anm can handle this many particles in a reasonable (e.g. days or a few weeks) run time [33].
Therefore, this section presents two examples, one of coarse aggregates in a mortar matrix, and one of sand in a cement paste matrix. A normal concrete mix is given in Table 1, in terms of mass used per cubic meter of concrete for both fine and coarse aggregates, which will be used to simulate mortar and concrete separately using the Anm model. The grading of the coarse and fine aggregates is given in the lower portion of Table 1.
In the simulations to be described next, the concrete specimen is a cube with side length 150 mm, filled with coarse aggregates, and the mortar specimen is in the shape of a cube with edge length 10 mm, filled with fine aggregates. The values of the many parameters of the Anm model were chosen to be able to create the desired structures. However, no attempt was made to optimize these parameter values-their values were set high in order to make sure the desired aggregate concentrations were achieved. The various internal parameters in the Anm model need to be systematically investigated as a function of sieve ranges used to see how these parameters can be optimized to achieve a wide range of particle-based microstructures. A first attempt at this systematic investigation will be carried out using the 2nd version of the Anm model, employing the more efficient contact function and partial parallelization to decrease run time [33].

Mortar example
Periodic material boundary conditions were employed for the mortar model. The total mass of sand used in the mortar specimen was 1.397 g, with volume of 527 mm 3 , and taking up 52.7 % of the specimen volume. The mass in each sieve range, 0.125 mm to 4 mm, is found by scaling the mass amounts in Table 1 by the factor 1.397 g/983 kg. The result of the Anm model is illustrated in Fig. 4 as a VRML image. The sand particles were Ottawa sand as used in the ASTM C109 mortar cube test. On a typical desktop PC, the time needed to build this simulated mortar was about 11 weeks, for 26,683 particles. This is a long run time, but there were several reasons for this. First, the Anm code was not optimized in any way-either through exploiting parallelism or in choosing the internal parameters that control how often different scenarios (e.g. position, orientation, size) were tried for each particle placed. These parameters (the maximum allowed attempts to pick, rotate, scale, and place a particle and the number of points chosen to fix the fourth unknown in the overlap algorithm) were set to high values since we wanted to get the necessary sand placement for this demonstration example and were not concerned about the run time. The PC processor used ran at 3.20 GHz and was manufactured in 2009. The actual run times on more modern, faster Windows machines, along with code optimization, a degree of parallelism in some of the internal Anm algorithms, and more judicious selection of the Anm parameters mentioned above are more than a factor of 20 less [33].
To get an idea of how run times vary with the volume fraction to be achieved, using the present version of Anm, we ran the same system with the same sand particle size distribution but only up to 30 %  Fig. 5 using a VRML image. The run time for this concrete simulation was about 2 weeks, for 5,640 particles, which was almost twice as long as the 30 % volume fraction mortar but only for about one-third the particles. We observed that building composites using non-periodic boundary conditions seem to take somewhat longer, for systems with an equal number of particles, than with periodic boundary conditions. This is probably due to the difficulty of placing particles near the boundary, which removes some of the unit cell volume from placement possibility as compared to when using periodic boundary conditions. This time difference is probably also due to the fact that the concrete gravel size distribution is narrower than the mortar sand. A narrower particle size distribution tends to come up more quickly against the monosize particle jamming limit, which can significantly slow down the particle placement process [4].

Potential model applications and improvements
The potential applications of the present version of the Anm model are many. With certain improvements, stated below, the number of potential applications is increased still further. First, the Anm model can give realistic mortar or concrete aggregate geometric structures. One can mesh these models with a finite element mesh and go on to calculate mechanical properties, including fracture, once appropriate mechanical properties have been assigned to the matrix and inclusions. A meshing scheme especially adapted to the Anm model has already been developed [19]. Discrete element models that attempt to simulate the flow of non-spherical powders suffer from two limitations [18,34]: (1) it is difficult to find realistic 3D particles that can serve as the starting point for the usual DEM method of modeling irregular particles, with clumps of spheres, and (2) it is difficult to generate an initial configuration of particles that can be relaxed into a mechanically stable starting point for the DEM simulation. Using the Anm model, one could easily generate large systems of randomly-shaped particles that can be used as geometrical input into discrete element software both for generating individual sphere clump models and for generating initial configurations of many particles.
Computational rheology models [35] also need a similar starting point, with 3-D particles dispersed in a matrix, upon which to exert a shear force and compute time-dependent particle positions. The Anm model can be easily used to generate such starting structures. Finally, the Anm model can be used to numerically generate various aggregate structures for which geometric statistics can be generated [36], such as radial distribution functions (how many particles are within a certain distance from a given aggregate particle) or the Fig. 5 The composite material structure of the concrete specimen with irregular shape crushed stones. The computational cell is shown as a translucent cube distance distribution between particle surfaces. These geometric statistics can then be compared to experimental data from X-ray CT to better understand the real aggregate structure in mortar and concrete, and the effect of real, non-spherical shape, which is a major variable in many mortar and concrete properties.
The interfacial transition zone exists between aggregate particles and the cement paste matrix [37,38]. This has been modeled as a uniform-thickness shell placed around each aggregate particle [39]. Placing a uniform-thickness shell around a spherical aggregate particle is simple, since it is just formed by a slightly larger spherical shell concentric with the original particle. This process is much harder for an irregular particle. Recently, a method has been developed that allows such a shell to be placed around a spherical harmonic particle by recalculating and slightly changing the spherical harmonic coefficients [32]. By combining this new algorithm with the Anm model, interfacial transition zones could be represented in these simulated composite structures, for the expansion of its ability to simulate real mortar and concrete [33].

Summary
In this paper, a geometrical model was described that simulates the geometrical structures of cementitious materials: the Anm model. Any arbitrary star-shaped particle can be represented in terms of spherical harmonic coefficients a nm , which will allow any sand or gravel used in construction to be employed in a simulation, giving the Anm algorithm broad applicability. Sand in mortar or gravel in concrete can be regarded as irregular shape particles in a matrix. The shapes of a class of aggregate can be captured by an X-ray CT scan, and then represented by spherical harmonic expansion, with the coefficients for each particle stored in a shape database. Such a database has been included in the Anm model, with five different kinds of particle in it, including sand, gravel, and portland cement, and a method has been demonstrated for generating new particles based on the shape statistics of real particles. The composite geometrical structure can be reproduced by parking multiple particles with the appropriate shape characteristics into an empty container. No overlap is allowed between aggregates, which is checked using a particle contact function. Both periodic and non-periodic material boundaries are implemented in the Anm model, and examples of both, for mortar and concrete, were described. Using Anm results as a computational input, the simulated composite structures of mortar and concrete can be further analyzed in terms of geometry, mechanical performance (e.g. elastic moduli, fracture paths), and transport properties (e.g. conductivity and diffusivity), and can be used as an initial configuration for DEM and rheology computations.