Acta Polytechnica Vol. 52 No. 6/2012 Geometrical Modeling of Concrete Microstructure for the Assessment of ITZ Percolation Daniel Rypl1, Tomáš Bým2 1Department of Mechanics, Faculty of Civil Engineering, Czech Technical University in Prague, Thákurova 7, 166 29 Prague, Czech Republic 2Golder Associates AB, Östgötagatan 12, 116 25 Stockholm, Sweden Corresponding author: drypl@fsv.cvut.cz Abstract 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 representing ITZ. Keywords: Percolation, concrete, interfacial transition zone, aggregates, hard core – soft shell, spherical harmonic analysis. 1 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 mi- crostructure. 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 mate- rials evolves in time, so that the percolation transition occurs after an ageing time (as it is typical for cement- based 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 ce- ment 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 com- putationally not simple, as the geometry and topology of this phase is complex. Note that micrometer and millimeter length scales have to be considered simul- taneously 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 uti- lized. In this continuum model, each randomly-placed 38 Acta Polytechnica Vol. 52 No. 6/2012 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 aggre- gate volume fraction was found to be dependent on the aspect ratio of the (ellipsoidal) aggregates. It is interesting, however, that the ITZ percolation thresh- old 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 three- dimensional 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 assess- ment 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. 2 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 direc- tion 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 x(η,ϕ) = r(η,ϕ) sin(η) cos(ϕ) y(η,ϕ) = r(η,ϕ) sin(η) sin(ϕ) (1) z(η,ϕ) = r(η,ϕ) cos(η) . η y x ϕ O z η ϕr ( , ) η y x ϕ O z η ϕr ( , ) surface surface real digital Figure 1: Description of the aggregate particle in the spherical coordinate system. 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 r(η,ϕ) = ∞∑ n=0 n∑ m=−n anmY m n (η,ϕ), (2) where anm are the coefficients of the expansion and Y mn (η,ϕ) are the spherical harmonic functions (see [14] for details). The coefficients anm are determined by the integration anm = ∫ 2π 0 ∫ π 0 r(η,ϕ) sin(η)Ȳ mn (η,ϕ) dη dϕ, (3) where Ȳ mn (η,ϕ) is the complex conjugate to Y mn (η,ϕ). An analytical evaluation of anm is practically not feasible (except for the case of a sphere where a00 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 r(η,ϕ) = N∑ n=0 n∑ m=−n anmY m n (η,ϕ), (4) where N is the order of the expansion, there is a reduc- tion in the order of the numerical integration needed for an accurate evaluation of coefficients anm and con- sequently 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 ac- quainted 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 contribu- tions of the expansion terms for n > 20 are usually 39 Acta Polytechnica Vol. 52 No. 6/2012 negligible, and that order 128 of Gaussian numerical integration is sufficient in most cases. It is also im- portant to realize that growth of the expansion order may lead not to expected improvement in the geomet- rical 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 repre- sented 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 expan- sion, by the order of the expansion, and by the set of expansion coefficients. For an evaluation of most geo- metrical 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. 3 Packing Algorithm In order to model a concrete specimen with the ag- gregate 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 aggre- gate particles into a representative volume. This kind of algorithm should (i) ensure sufficient randomness, (ii) comply with the volumetric content and the sta- tistical distribution of size grades of the individual aggregates, and (iii) prevent overlapping of individ- ual 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 par- ticular 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 repre- sentative 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 1With respect to the present work only the volume, needed for the evaluation of aggregates volume fraction, and the first gradients of the r with respect to η and ϕ, needed for the intersection check, are of interest. 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). 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 ran- domly 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 bound- ary 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 grad- ing, alternative more sophisticated packing procedures (see [24, 25, 26] for example) are necessary. The crucial aspect of the packing algorithm is the in- tersection check. In order to make it efficient, an octree data structure is used to easily identify pairs of parti- cles which could potentially overlap. Although many intersection checks can be resolved using the bound- ing 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 estab- lished. The spheres in the refined bounding box are spheres circumscribed to individual simplices in a con- 40 Acta Polytechnica Vol. 52 No. 6/2012 a) b) c) 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. 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, how- ever, the basic triangulation may be too coarse and must be globally refined by a regular recursive subdi- vision, in which each surface triangle is divided into four geometrically similar subtriangles. The newly- introduced 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 in- dividual 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 single- level refinement, the refined bounding box comprises Figure 4: 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. 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 par- ticle 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 con- strained 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 expansion2) several al- most identical spheres (having a similar radius and location). Of course, this increases the size of the re- fined 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 main- taining 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 2Note that expansion is accompanied by an appropriate shift of the center of the sphere to minimize necessary growth of its radius. The center, however, has to be kept safely inside the particle. 41 Acta Polytechnica Vol. 52 No. 6/2012 a) b) c) Figure 5: Surface triangulation of an aggregate par- ticle: a) basic triangulation (240 triangles), b) tri- angulation with one-level refinement (960 triangles) and c) triangulation with two-level refinement (3840 triangles). a) b) c) d) e) f) Figure 6: Refined bounding box (of the particle from Figure 5): a) no merging applied (240 spheres), b) merging using expansion factor 0.05 (156 spheres), c) 0.1 (71 spheres), d) 0.15 (49 spheres), e) 0.2 (36 spheres) and f) 0.25 (31 spheres). ray length r(η,ϕ), the refined bounding box may be precomputed for individual particles in the container and then only appropriately scaled and rotated ac- cording to the randomly generated values during the packing procedure. The refined bounding box is used to efficiently resolve the vast majority of the inter- section 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 perpen- dicular 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 EC 1 90 2min d(MS , S ) 2min d(MS , S ) EC 2 1 1 1 MS1 2 1 2 MS 12 2 d 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. ECi indicates expansion center of particle i, MSji denotes marginal sphere (in dark gray) of particle i facing particle j, d(MSji ,Sj) 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(MS21,S2)) and min(d(MS12,S1)). (on the side facing the other particle) shifted toward the other particle by the so far detected minimum dis- tance 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 stag- gered projection may become quite inefficient if the 3Note that the refined bounding box is generally not suitable to confirm the intersection of two particles because the thickness of the cover of the particle by the refined bounding box is variable. 42 Acta Polytechnica Vol. 52 No. 6/2012 EC 2 2 1EC 1 SP A A B B1 1 2 2 SPA B Figure 8: 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 SPA as starting point results in locally clos- est points A1 and A2. Globally closest points B1 and B2 are obtained by choosing SPB as starting point. surfaces of the particles facing each other are almost parallel, or if the particles are far from each other. 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 project- ing 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 A1 and A2 in Figure 8), depending on the starting point. In order to identify the closest points in the global sense (points B1 and B2 in Figure 8), the staggered projec- tion 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 refine- ment 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 4There is a pathological case (not too likely to happen because of round-off errors), when the starting point is located on a line which is normal simultaneously to both overlapping particles. In this case, the closest points are distinct, each one inside the other particle. 1EC 1 EC 2 2 1 2 BS BS Figure 9: Points (white circles) used as starting points for the intersection check of particles repre- sented by the real geometry. BSi indicates the single bounding sphere of the i-th particle. 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 num- ber 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 scal- ing 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 al- gorithm, 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 surround- ing 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 esti- mating 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 43 Acta Polytechnica Vol. 52 No. 6/2012 a)b) d)c) Figure 10: Aggregates packing procedure: a) con- tainer of aggregate shapes (particles are rotated accord- ing to the final location; only used shapes are shown), b) empty representative volume, c) representative vol- ume with periodic aggregates pack, d) representative volume with periodic expanded aggregates pack. 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, de- pending on the order in which the aggregate particles are processed, and whether a maximum scaling factor (cumulative within a single expansion step) is pre- scribed. Currently, the simplest approach is adopted, in which the aggregate particles are processed in the order in which they were packed into the representa- tive 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 two- dimensional example of a periodic representative vol- ume 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. 4 Assessment of ITZ Percolation After the representative aggregate pack is available, the ITZ percolation has to be verified. This is accom- Figure 11: Some of the representative shapes of aggregate particles. 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 aggre- gates 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 intercon- nects 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 correspond- ing 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. 5 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 44 Acta Polytechnica Vol. 52 No. 6/2012 Figure 12: Aggregates packing into representative cube (levels of the gray correspond to individual size grades). 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 inte- gration. 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. 6 Conclusions This paper has presented an algorithm for geometrical modeling of the microstructure of concrete with em- bedded aggregate particles of realistic shape. Unlike Figure 13: Detail of aggregates packing into repre- sentative cube. 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 ran- domly 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 alterna- tive packing procedures that enable higher aggregate volumetric content values to be reached. Acknowledgements This work was supported by the Ministry of Educa- tion of the Czech Republic under project No. MSM 6840770003. Its financial assistance is gratefully ac- knowledged. 45 Acta Polytechnica Vol. 52 No. 6/2012 References [1] Garboczi, E.J., Bentz, D.P.: Computer Simula- tion and Percolation Theory Applied to Concrete, Annual Reviews Of Computational Physics VII, D. Stauffer (ed.), World Scientific Publishing Company, 1999, pp. 85–123. [2] Bentz, D.P., Garboczi, E.J.: Percolation of Phases in a Three-Dimensional Cement Paste Microstruc- ture Model, Cement and Concrete Research, 21, 1991, pp. 325–344. [3] Bentz, D.P., Schlangen, E., Garboczi, E.J.: Com- puter Simulation of Interfacial Zone Microstruc- ture and its Effect on the Properties of Cement- Based Composites, in Material Science of Con- crete IV, Skalny, J.P. et al. (eds), Journal of the American Ceramic Society, 1995, pp. 155–199. [4] Garboczi, E.J., Bentz, D.P.: Computer Simula- tion of the Diffusivity of Cement-Based Materials, Journal of Materials Science, 27, 1992, pp. 2083– 2092. [5] Shane, J.D., Mason, T.O., Jennings, H.M., Gar- boczi, E.J., Bentz, D.P.: Effect of the Interfacial Transition Zone on the Conductivity of Portland Cement Mortars, Journal of the American Ce- ramic Society, 83, 2000, pp. 1137–1144. [6] Maso, J.C.: The Bond Between Aggregates and Hydrated Cement Paste, in Proceedings of the 7th International Congress on the Chemistry of Cement I, VII-1/3, 1980. [7] Garboczi, E.J., Bentz, D.P.: Digital Simulation of the Aggregate-Cement Paste Interfacial Zone in Concrete, Journal of Materials Research, 6, 1991, pp. 196–201. [8] Garboczi, E.J.: Finite Element and Finite Differ- ence Programs for Computing the Linear Electric and Elastic Properties of Digital Images of Ran- dom Materials, NIST Internal Report 6269, 1998. [9] Jia, X., Williams, R.A.: A Packing Algorithm for Particles of Arbitrary Shapes, Powder Technology, 120, 2001, pp. 175–186. [10] Tourquato, S.: Two-point Distribution Function for a Dispersion of Impenetrable Spheres in a Ma- trix, Journal of Chemical Physics, 85, 1986, pp. 6248–6249. [11] Lee, S.B., Torquato, S.: Porosity for the Penetrable-Concentric-Shell Model of Two-Phase Disordered Media: Computer Simulation Results, Journal of Chemical Physics, 89, 1988., pp. 3258– 3263 [12] Bentz, D.P., Garboczi, E.J., Stutzman, P.E.: Computer Modelling of the Interfacial Transition Zone in Concrete, in Interfaces in Cementitious Composites, Maso, J.C. et al. (eds), 1993, pp. 259–268. [13] Bentz, D.P., Hwang, J.T.G., Hagwood, C., Garboczi, E.J., Snyder, K.A., Buenfeld, N., Scrivener, K.L.: Interfacial Zone Percolation in Concrete: Effects of Interfacial Thickness and Aggregate Shape, in Microstructure of Cement- Based System/Bonding and Interfaces in Cemen- titious Materials, Diamond, S. et al. (eds), Ma- terial Research Society, Pittsburgh, 1995, pp. 437–442. [14] Garboczi, E.J.: Three-dimensional Mathematical Analysis of Particle Shape Using X-ray Tomog- raphy and Spherical Harmonics: Application to Aggregate Used in Concrete, Cement and Con- crete Research, 32, 2002, pp. 1621–1638. [15] Berryman, J.G.: Random Close Packing of Hard Spheres and Disks, Physical Review A, 27, 1983, pp. 1053–1061. [16] Salvat, W.I., Mariani, N.J., Barreto, G.F., Martínez, O.M.: An Algorithm to Simulate Pack- ing Structure in Cylindrical Containers, Catalysis Today, 107–108, 2005, pp. 513–519. [17] Kristiansen, K.L., Wouterse, A., Philipse, A.: Simulation of Random Packing of Binary Sphere Mixtures by Mechanical Contraction, Physica A, 358, 2005, pp. 249–261. [18] Knott, G.M., Jackson, T.L., Buckmaster, J.: The Random Packing of Heterogeneous Propellants, AIAA Journal, 39, 4, 2001, pp. 678–686. [19] Lo, S.H., Wang, W.X.: Generation of Tetrahedral Mesh of Variable Element Size by Sphere Pack- ing over an Unbounded 3D Domain, Computer Methods in Applied Mechanics and Engineering, 194, 2005, pp. 5005–5018. [20] Han, K., Feng, Y.T., Owen, D.R.J.: Sphere Pack- ing with a Geometric Based Compression Al- gorithm, Powder Technology, 155, 1, 2005, pp. 33–41. [21] Bažant, Z.P., Tabbara, M.R., Kazemi, M.T., Pijaudier-Cabot, G.: Random Particle Model for Fracture of Aggregate of Fiber Composite, Journal of Engineering Mechanics - ASCE, 116, 8, 1990, pp. 1686–1705. [22] Wang, Z.M., Kwan, A.H.K., Chan, H.C.: Meso- scopic Study of Concrete I: Generation of Random Aggregate Structure and Finite Element Mesh, Computers and Structures, 70, 5, 1999, pp. 533– 544. 46 Acta Polytechnica Vol. 52 No. 6/2012 [23] Koenke, C., Eckardt, S., Haefner, S.: Spatial and Temporal Multiscale Simulations of Damage Processes for Concrete, in Innovation in Compu- tational Structures Technology, Topping B.H.V. et al. (eds), Saxe-Coburg Publications, 2006, pp. 133–157. [24] Van Mier, J.G.M., van Vliet, M.R.A.: Influence of Microstructure of Concrete on Size/Scale Ef- fects in Tensile Fracture, Engineering Fracture Mechanics, 70, 16, 2003, pp. 2281–2306. [25] Leite, L.P.B., Slowik, V., Mihashi, H.: Computer Simulations of Fracture Processes of Concrete using Meso-Level Models of Lattice Structures, Cement and Concrete Research, 34, 2004, pp. 1025–1033. [26] Haefner, S., Eckardt, S., Koenke, C.: A Geo- metrical Inclusion-Matrix Model for the Finite Element Analysis of Concrete at Multiple Scales, in Proceedings of the 16th IKM, Guerlebeck, K. et al. (eds), 2003. [27] Baker, T.J.: Automatic Mesh Generation for Complex Three-dimensional Regions Using a Con- strained Delaunay Triangulation, Engineering with Computers, 5, 1989, pp. 161–175. [28] Bažant, Z.P., Caner, F.C., Carol, I., Adley, M.D., Akers, S.A.: Microplane Model M4 for Concrete. Part I: Formulation with Work-Conjugate Devia- toric Stress, Journal of Engineering Mechanics - ASCE, 126, 2000, pp. 944–953. [29] Rypl, D.: Discretization of Three-Dimensional Aggregate Particles, in CD-ROM Proceedings of the III. European Conference on Computational Mechanics, Mota Soares, C.A. et al. (eds), 2006. 47