Geometrical Modeling of Concrete Microstructure for the Assessment of ITZ Percolation

Percolation is considered to be a critical factor affecting the transport properties of multiphase materials. In the case of concrete, the transport properties are strongly dependent on the interfacial transition zone (ITZ), which is a thin layer of cement paste next to aggregate particles. It is not computationally simple to assess ITZ percolation in concrete, as the geometry and topology of this phase is complex. While there are many advanced models that analyze the behavior of concrete, they are mostly based on the use of spherical or ellipsoidal shapes for the geometry of the aggregate inclusions. These simplified shapes may become unsatisfactory in many simulations, including the assessment of ITZ percolation. This paper deals with geometrical modeling of the concrete microstructure using realistic shapes of aggregate particles, the geometry of which is represented in terms of spherical harmonic expansion. The percolation is assessed using the hard core – soft shell model, in which each randomly-placed aggregate particle is surrounded by a shell of constant thickness


Introduction
Concrete is nowadays a modern composite material.It is a multiscale material with length scales from nanometers (C-S-H), via micrometers (cement paste) to millimeters (mortar and concrete).It is a different random composite at each length scale.Thus the range of the microstructure of concrete spans nine orders of magnitude.It is therefore a large and difficult task to try to relate the microstructure and the properties of concrete theoretically.
A possible way to relate several material properties to the microstructure of the material are based on percolation theory [1,2].Percolation theory is using the idea of connectivity.The percolation threshold denotes the volume fraction of a particular phase of a composite material at which that phase goes from being disconnected to being connected (or vice versa), so that there is a change in the topology of the microstructure.Percolation properties are now more or less commonly accepted as the critical geometrical and topological factors influencing the transport properties of multiphase materials (e.g.ionic diffusivity, electric or thermal conductivity) [3,4,5].Moreover, in many practical applications the structure of composite materials evolves in time, so that the percolation transition occurs after an ageing time (as it is typical for cementbased materials and gels).These are very important factors, because it is now well recognized that much of the deterioration of concrete in infrastructures is caused by corrosion of the reinforcing bars due to a massive chloride attack (or due to an attack by ions of some other salts).
In the case of mortar and concrete, the transport properties are strongly dependent on the region of cement paste close to the aggregate particle surface (typically within 50 micrometers).This region, known as the interfacial transition zone (ITZ), exhibits higher capillary porosity and larger pores than the bulk cement paste matrix [6,7].These features are commonly attributed to the cement particle packing effect and the one-side growth effect.However, if ITZs do not percolate, their effect on transport will be fairly small, as any transport path through concrete would have to go through the bulk cement paste.The transport properties would then be dominated by the transport properties of the bulk cement paste.
The problem of ITZ percolation in concrete is computationally not simple, as the geometry and topology of this phase is complex.Note that micrometer and millimeter length scales have to be considered simultaneously in such a study.This makes the standard models based solely on digital image processing [8,9] prohibitive in terms of memory requirements.The hard core -soft shell model [10,11] is therefore utilized.In this continuum model, each randomly-placed aggregate particle is surrounded by a shell of constant thickness representing the ITZ.While the hard core aggregates may not overlap one another, the soft shell ITZs are free to overlap one another.ITZ percolation is verified if it propagates continuously through the whole specimen in a particular direction.Until now, the ITZ percolation problem has been studied for spherical [12] and ellipsoidal aggregate particles [13].The ITZ percolation threshold in terms of the aggregate volume fraction was found to be dependent on the aspect ratio of the (ellipsoidal) aggregates.It is interesting, however, that the ITZ percolation threshold in terms of the volume fraction of the ITZs was not dependent on the aggregate aspect ratio [13].
This work goes one step further and elaborates the assessment of ITZ percolation for realistic shapes and distribution of aggregates.A relatively simply and robust approach introduced in [14] is employed for describing real three-dimensional aggregate particles.This method is based on approximating the particle shape by spherical harmonic functions (the threedimensional equivalent of two-dimensional Fourier analysis).Although this representation is not universal (it implies that the aggregate particle is star-like in shape with no internal voids), it is suitable for almost all aggregates used in structural concrete.A significant advantage of this approach is that the resolution of the smooth representation can be flexibly controlled by the number of terms in the expansion.
The paper is organized as follows.In Section 2, the representation of the aggregate particle using spherical harmonic expansion is recalled.Then the algorithm for packing the aggregate particles into the representative volume is described in Section 3. The actual assessment of ITZ percolation is worked out in Section 4. Section 5 presents a simple example of a concrete specimen with realistic shapes and distribution of the aggregates.The paper ends with concluding remarks in Section 6.

Geometrical Representation of Aggregate Particles
The geometrical representation of aggregate particles is based on a scalar function r(η, ϕ), defined as the distance of the particle surface point from the particle center (possibly center of mass) measured in the direction of spherical coordinates η and ϕ with the origin located at that center (Figure 1).The Cartesian coordinates of the surface point (components of its position vector r) are then given by If the function r(η, ϕ) is smooth on unit sphere (0 ≤ η ≤ π, 0 ≤ ϕ ≤ 2π) and periodic in ϕ, it can be expressed in the form of the expansion into spherical harmonic functions as where a nm are the coefficients of the expansion and Y m n (η, ϕ) are the spherical harmonic functions (see [14] for details).The coefficients a nm are determined by the integration where Ȳ m n (η, ϕ) is the complex conjugate to Y m n (η, ϕ).An analytical evaluation of a nm is practically not feasible (except for the case of a sphere where a 00 is the only non-zero coefficient), and it is therefore necessary to apply numerical integration.Gaussian numerical integration is employed in the present study.Theoretically, the order of the numerical integration should approach infinity.However, taking for the approximation of r(η, ϕ) the final summation where N is the order of the expansion, there is a reduction in the order of the numerical integration needed for an accurate evaluation of coefficients a nm and consequently also in the number of values r(η, ϕ) needed for the integration.The values of r(η, ϕ) at the integration points may be derived, for example, from a digital (voxel-based) representation of an aggregate particle (easily acquainted from computer tomography or any other similar scanning device).r(η, ϕ) is defined as the length of the segment in the direction given by η and ϕ connecting the particle center with the side of the voxel forming the aggregate particle boundary (Figure 1 on the right).
Numerical experiments reveal that the contributions of the expansion terms for n > 20 are usually negligible, and that order 128 of Gaussian numerical integration is sufficient in most cases.It is also important to realize that growth of the expansion order may lead not to expected improvement in the geometrical representation, but to undesirable capturing of the unrealistic digital roughness inherently present in the voxel-based description.An example of the representation of a particular aggregate particle using spherical harmonic expansion with 128-point Gaussian numerical integration is shown in Figure 2 for different values of N .
After the spherical harmonic analysis is completed, the aggregate particle in a general location is represented by its origin and unit vectors, corresponding to the positive x, y, and z axes of the local spherical coordinate system used for spherical harmonic expansion, by the order of the expansion, and by the set of expansion coefficients.For an evaluation of most geometrical properties1 only the last two parameters, the order of the expansion and the expansion coefficients, are relevant.The remaining parameters describe the spatial location of the aggregate particle, and may be used for evaluating some of the geometrical properties in the global Cartesian coordinate system.

Packing Algorithm
In order to model a concrete specimen with the aggregate distribution reasonably matching the realistic shapes and volumetric fractions of aggregate particles of various size grades, it is necessary to employ an appropriate packing algorithm that assembles aggregate particles into a representative volume.This kind of algorithm should (i) ensure sufficient randomness, (ii) comply with the volumetric content and the statistical distribution of size grades of the individual aggregates, and (iii) prevent overlapping of individual aggregate particles.This is a non-trivial task when realistic aggregate shapes are considered.There are many packing strategies [15,16,17,18,19,20] designated for spherical, cylindrical, and ellipsoidal inclusions, but many of them are specialized for a particular purpose, or are restricted by some constraints not applicable to the problem of aggregate packing.The most common packing approach is based on the so-called take-and-place method ( [21,22,23]), which is also employed in the present study.
In the first phase, the particles are randomly chosen from a container of predefined and sufficiently representative shapes of aggregates of a particular type (e.g.artificially crushed aggregates or natural aggregates from a riverbed).Each particle is randomly scaled to fit into a particular size grade and, its volume is added to the volumetric content of that grade.The individual size grades are processed sequentially in descending order with respect to the sieve size, until the desired volumetric content is reached.In the next phase, which is computationally most demanding, the individual particles are placed into the representative volume in descending order with respect to the radius of their bounding sphere.Each particle is first randomly rotated, and then its tentative location in the representative volume is randomly generated.Note that if the periodic representative volume is treated, the periodic counterpart(s) of the processed particle must also be handled.Before the tentative location is accepted, it has to be verified that the particle (and its periodic counterpart(s)) does not overlap with any of the already inserted particles and, if the periodicity is not considered, that it does not intersect the boundary of the representative volume.If an overlap or an intersection is detected, a new position (and after a large number of unsuccessful attempts also a new rotation) is randomly generated and the process is repeated until the location can be accepted.Note that the above algorithm is suitable only for realistic size grading with volumetric content of the aggregates not exceeding 50-60 % of the total representative volume.For higher volumetric content or for unusual size grading, alternative more sophisticated packing procedures (see [24,25,26] for example) are necessary.
The crucial aspect of the packing algorithm is the intersection check.In order to make it efficient, an octree data structure is used to easily identify pairs of particles which could potentially overlap.Although many intersection checks can be resolved using the bounding spheres of individual particles (a single bounding sphere per particle), the computational effort related to the remaining checks is still prohibitive.A refined bounding box, comprising a set of spheres of variable radius covering the whole particle, is therefore established.The spheres in the refined bounding box are spheres circumscribed to individual simplices in a con- strained 3D Delaunay triangulation [27] constructed over a set of particle surface points.These points are obtained as the intersection of the particle surface with regularly spaced rays emanating from the expansion center of the particle.The particle surface points are connected by a surface triangulation, the topology of which is the same for all particles and is known in advance.The process for constructing the refined bounding box is schematically demonstrated for a 2D case in Figure 3.In the current implementation, the directions of the rays are taken from the microplane material model [28], which uses 122 regularly spaced directions (defining normals of individual constitutive microplanes in that material model).This yields 240 triangles approximating the particle surface.This basic triangulation is sufficient up to expansion order 10 (inclusive).For higher expansion order values, however, the basic triangulation may be too coarse and must be globally refined by a regular recursive subdivision, in which each surface triangle is divided into four geometrically similar subtriangles.The newlyintroduced subtriangle vertices define new rays which, in turn, produce new surface points that enrich the original set of basic surface points, and that replace the original vertices (obtained by subdivision) of individual subtriangles (see Figure 4).An example of particle surface triangulations corresponding to basic points and to points after the first and the second refinement level is depicted in Figure 5.The estimate of the number of refinement levels ensuring that the particle is completely covered by the refined bounding box is based on a heuristics (verified experimentally on a set of particles) which applies a one-level refinement (an increase in the number of triangles by a factor of 4) whenever the expansion order increases by 10.This implies that just one level of refinement is necessary for the most commonly-used expansion orders between 10 and 20 (inclusive).However, even for just a singlelevel refinement, the refined bounding box comprises 960 spheres, which is computationally unfavourable.
The number of spheres forming the refined bounding box can be decreased by performing the subdivision locally (for example, with respect to the change in particle surface curvature).In this case, the constrained Delaunay triangulation is no longer conforming, due to the introduced hanging nodes (if only one out of two adjacent triangles was refined).However, this causes no difficulty.The only consequence is that the refined bounding box is slightly larger because it is formed by fewer spheres.In the current approach, an alternative strategy is adopted to keep the number of bounding spheres reasonably small.Initially, the constrained Delaunay triangulation is built using only the basic 240 surface triangles, irrespective of the actual number of refinement levels that are applied.Then, for each of the created tetrahedra (a maximum of 240), its circumscribed sphere is appropriately expanded2 to cover also the refined points constructed during the refinement over the basic triangle(s) bounding that tetrahedron.Finally, the total number of spheres in the refined bounding box is reduced by merging (again using the appropriate expansion 2 ) several almost identical spheres (having a similar radius and location).Of course, this increases the size of the refined bounding box.The magnitude of the increase is controlled by an expansion factor related to the radius of the single bounding sphere of the particle.Figure 6 demonstrates the effect of merging the spheres in the refined bounding box of the particle from Figure 5 when one refinement level (see Figure 5b) was used.
In the present study, an expansion factor of 0.1 was chosen as optimal, considerably reducing the number of spheres in the refined bounding box while maintaining reasonably tight bounding.Unfortunately, it is difficult to make a quantitative assessment of the measure of agreement between the particle itself and its refined bounding box.
Since the particle is geometrically described by the ray length r(η, ϕ), the refined bounding box may be precomputed for individual particles in the container and then only appropriately scaled and rotated according to the randomly generated values during the packing procedure.The refined bounding box is used to efficiently resolve the vast majority of the intersection checks between close particles.Note that handling all pairs of spheres (one from each refined bounding box) in such a check is inefficient, even for the relatively small number of spheres in refined bounding boxes.In the present implementation, the most marginal sphere (on the side facing the other particle along the direction connecting the expansion centers of the two particles) in the refined bounding box of each of the two particles is identified first.Then a search is made for the closest sphere from the refined bounding box of one particle to the marginal sphere of the refined bounding box of the other particle.From these two pairs, the one with smaller distance between the spheres is selected.Then the search continues only among mutual pairs of those spheres (one from each refined bounding box) crossing a layer perpendicular to the line connecting the expansion centers of the two particles.The planes bounding the layer are defined as the planes touching the marginal spheres (on the side facing the other particle) shifted toward the other particle by the so far detected minimum distance between spheres defining the refined bounding boxes.The strategy discussed above is schematically depicted in Figure 7.If the intersection cannot be reliably refuted (i.e. if any of the spheres in the refined bounding box of one particle is touching or overlapping any sphere in the refined bounding box of the other particle) or cannot be reliably proved3 (if the center of any sphere in the refined bounding box of one particle is inside the other particle), then the real geometry of the investigated particles is considered.Since the surface of the particle (represented by the spherical harmonic expansion) is naturally parametrized by the spherical angles η, ranging from 0 to π and ϕ, varying from 0 to 2π and being periodic, a standard algorithm (the so-called closest point projection) is adopted which returns, for a given point, the closest point on the surface (for details see [29]).To obtain a pair of mutually closest points, one on each particle surface, the closest point projection is performed in a staggered way.This means that the projection of a point to the surface of one particle is then projected to the surface of the other particle, and so on (see Figure 8), until the converged state is achieved (the projection of the point on the surface of one particle to the surface of the other particle is the point on the other particle, and vice versa).Note that the staggered projection may become quite inefficient if the EC 2 In the case of convex non-overlapping particles, this procedure yields the closest points in the global sense, irrespective of the location of the starting point.If the convex particles overlap or touch each other, then the procedure converges to geometrically identical points on the intersection of the particles4 .To speed up identification of the overlapping particles, the points projected to the surface of one particle are checked against being inside the other particle before projecting it to the other particle.If the particles are not convex (which is the case for the present study), the staggered projection generally results in points closest to each other in the local sense only (points A 1 and A 2 in Figure 8), depending on the starting point.In order to identify the closest points in the global sense (points B 1 and B 2 in Figure 8), the staggered projection procedure is performed for several starting points.
In the present implementation, the points used to build the refined bounding box are used.To limit the total number of starting points, only those points on the surface of one particle which are inside the single bounding sphere of the other particle are considered (see Figure 9).From this point of view, local refinement of the basic surface triangulation (used to build the refined bounding box) has a clear advantage over global refinement, as it yields fewer points that can potentially be used as starting points for the staggered projection.The performance can also be improved by using only those surface points of one particle, within the single bounding sphere of the other particle, which are closest to any of the surface points on the other particle.However, this acceleration may determine only the locally closest points.Note that, in some pathological cases, there may be no starting points or only a few starting points, within the single bounding spheres, even if the refined bounding boxes overlap.In this case, the radii of the single bounding spheres are slightly enlarged to identify a large enough number of starting points.It is apparent that the above approach is appropriate for identifying whether two particles do or do not intersect, and for determining the distance that is needed for estimating the scaling factor (see below) of close enough particles (with overlapping refined bounding boxes).The distance of two particles with non-overlapping refined bounding boxes is only approximated by the minimum distance of the spheres in the refined bounding boxes.If their distance is to be evaluated exactly, then the starting points for the staggered projection are determined using enlarged single bounding spheres containing the expansion center of the other particle.
In the second phase of the aggregate packing algorithm, each of the aggregate particles is expanded (scaled with respect to its expansion center), and is shifted until it is touching at least two of its surrounding particles or, if periodicity is not considered, one or more sides of the representative volume.Note that this may lead to the introduction of additional periodic particles, if a representative volume with a periodic boundary is considered.The expansion is a computationally demanding process that must be performed in an iterative manner.In order to reduce the computational load related to the calculation of the distance of two aggregate particles (needed for estimating the scaling factor) and also to the intersection check (needed for verifying that the scaling has not been overestimated), only particles in the immediate neighbourhood are considered.Since a large number of iterations are required to achieve touching (within the round-off error) of two neighbouring particles, the scaling yielding non-overlapping particles closer than the sum of the thicknesses of their ITZs is accepted as approximate touching.Note that although this approach allows us to handle variable ITZ thickness, we utilize a constant thickness that is the same for all particles.There are several expansion scenarios, depending on the order in which the aggregate particles are processed, and whether a maximum scaling factor (cumulative within a single expansion step) is prescribed.Currently, the simplest approach is adopted, in which the aggregate particles are processed in the order in which they were packed into the representative volume.Note that the expansion phase results in a slight violation of the initially prescribed size grading.
The whole process of the packing procedure is schematically depicted in Figure 10 for a twodimensional example of a periodic representative volume with aggregate particles corresponding to three size grades (distinguished by different levels of gray).Note that some of the particles appear in the pack more than once.These are the periodic particles, and they can easily be identified in Figure 10 as those crossing the boundary of the representative volume.

Assessment of ITZ Percolation
After the representative aggregate pack is available, the ITZ percolation has to be verified.This is accom- plished using the hard core -soft shell model, in which the hard core is formed by the aggregate particle and the soft shell by ITZ of constant thickness (the same for all particles, irrespective of their size).Initially the connectivity of aggregate particles that are less than twice the thickness of ITZ apart is built up.Note that this connectivity has already been established during the aggregates expansion process of the aggregates packing procedure described in Section 3. The mutually connected aggregate particles form regions that are surrounded by distinct continuous ITZs.ITZ percolation occurs if any of the distinct ITZs interconnects the opposite sides of the representative volume (preferably in all principal directions).To prove that a particular pack has percolated, it is sufficient to tag aggregate particles touching (and/or crossing, if the periodicity of the pack is taken into account) opposite sides of the representative volume (using a different tag for each side) and to verify whether in the connectivity lists, there are particles marked by tags corresponding to opposite sides.Since it is difficult to evaluate the ITZ volume fraction, as ITZ is free to overlap itself, ITZ percolation is indirectly described by the aggregate volume fraction of the total representative volume at which ITZ percolation just occurs.Note that generally only the volume of the parts of the aggregate particles that are inside the representative volume should be taken into consideration.However, if the aggregate pack is built as periodic, then the whole volume of each particle is taken, with the exception of the periodic counterparts, which are ignored.

Example
In this section, a single aggregate pack used for ITZ percolation assessment is presented.The pack was generated inside a representative cube 50 mm in edge length, with a total volumetric content of aggregates prescribed to 50 %.A total of three aggregates size grades were used.The coarsest grade ranges from 8 to 16 mm, the intermediate grade from 4 to 8 mm, and the finest grade from 0.5 to 4 mm.The thickness of ITZ was chosen to be 30 µm.The container of representative aggregate shapes was formed by 100 aggregate particles derived from a randomly generated ellipsoid by randomly scaling the length of rays r used for evaluating the coefficients of the spherical harmonic expansion (Eq.( 3)).Note however that sufficient correlation was enforced for length of rays corresponding to adjacent integration points, in order to avoid unrealistic artificial shapes.The expansion coefficients were calculated for an expansion order equal to 20, using 128-point numerical Gaussian integration.Some of the representative aggregate particles are visualized in Figure 11.Although the particles are generated artificially, they quite realistically resemble natural riverbed aggregates.The final aggregate pack, shown in Figures 12 and 13, contains 2542 particles, out of which 14 are from the coarse grade, 161 from the intermediate grade, and 2367 are from the fine grade.The resulting aggregate volumetric content was 53 %, which suggests that there was not much space for particle expansion.However, the percolation of ITZ in this particular case was not proved.This again reveals that a more powerful packing procedure is desirable for the higher aggregate volumetric content at which ITZ percolation may occur.

Conclusions
This paper has presented an algorithm for geometrical modeling of the microstructure of concrete with embedded aggregate particles of realistic shape.Unlike other works on similar topics, real-shaped aggregate particles, with the geometrical representation based on spherical harmonic expansion have been considered.The microstructure that was produced is used for assessing ITZ percolation within the representative volume, using the hard core -soft shell model.The representative volume is built using a simple two-phase randomized particle packing procedure.In the first phase, the particles are taken from the container of representative aggregate particle shapes, and are randomly placed into the representative volume.In the second phase, the particles are scaled and shifted until they come into mutual contact (within the tolerance given by thickness of ITZ).To make the investigation of the particle contacts more efficient, the individual particles are approximated by a refined bounding box composed of a set of spheres of variable radius.If the contact cannot be reliably assessed using the refined bounding box, the real geometry of the investigated particles is considered.The contacts are used to set up the connectivity of the aggregate particles, which is then used to verify whether ITZ percolation occurs.In order to estimate the ITZ percolation threshold, given in terms of the aggregate volume fraction of the total representative volume, a large number of simulations have to be performed.Those that are close to the state at which ITZ percolation just occurs are then taken into consideration.However, this is a subject for future work.Future research will also focus on alternative packing procedures that enable higher aggregate volumetric content values to be reached.

Figure 1 :
Figure 1: Description of the aggregate particle in the spherical coordinate system.

Figure 2 :
Figure 2: Representation of an aggregate particle by spherical harmonic expansion of order 5, 10, 15, 20, 25, and 30 (from left to right and from top to bottom).

Figure 3 :
Figure 3: Construction of refined bounding box of an aggregate particle: a) surface points generated by the intersection of regularly spaced rays emanating from the expansion center of the particle with the surface of the particle, b) constrained Delaunay triangulation of surface points, c) refined bounding box defined by the envelope (in gray) of circles circumscribed to Delaunay simplices.

Figure 4 :
Figure4: Refinement of the surface triangulation using hierarchical subdivision.Black nodes correspond to surface points on the preceding refinement level, white nodes are subdivision points (not on the real surface) on current refinement level, gray nodes are surface points on the current refinement level.

Figure 7 :
Figure 7: Intersection check of refined bounding boxes of two particles.Only pairs of spheres (one from each refined bounding box) crossing the light gray layer are investigated.EC i indicates expansion center of particle i, M S ji denotes marginal sphere (in dark gray) of particle i facing particle j, d(M S j i , S j ) represents distance between marginal sphere of particle i and any sphere of particle j, d stands for the minimal detected distance (monotonically decreasing during the intersection check) initially set to smaller from min(d(M S2  1 , S 2 )) and min(d(M S 1 2 , S 1 )).

Figure 9 :
Figure 9: Points (white circles) used as starting points for the intersection check of particles represented by the real geometry.BS i indicates the single bounding sphere of the i-th particle.

Figure 10 :
Figure 10: Aggregates packing procedure: a) container of aggregate shapes (particles are rotated according to the final location; only used shapes are shown), b) empty representative volume, c) representative volume with periodic aggregates pack, d) representative volume with periodic expanded aggregates pack.

Figure 11 :
Figure 11: Some of the representative shapes of aggregate particles.

Figure 12 :
Figure 12: Aggregates packing into representative cube (levels of the gray correspond to individual size grades).

Figure 13 :
Figure 13: Detail of aggregates packing into representative cube.
Intersection check of particles represented by the real geometry using the closest point projection in the staggered way.Starting points are indicated as white circles, the intermediate projection points as gray circles, and the closest points as black circles.Choosing SP A as starting point results in locally closest points A 1 and A 2 .Globally closest points B 1 and B 2 are obtained by choosing SP B as starting point.
surfaces of the particles facing each other are almost parallel, or if the particles are far from each other.