3D morphological modeling of concrete using multiscale Poisson polyhedra

This paper aims at developing a random morphological model for concrete microstructures. A 3D image of concrete is obtained by micro-tomography and is used in conjunction with the concrete formulation to build and validate the model through morphological measurements. The morphological model is made up of two phases, corresponding to the matrix, or cement paste and to the aggregates. The set of aggregates in the sample is modeled as a combination of Poisson polyhedra of diﬀerent scales. An algorithm is introduced to generate polyhedra packings in the continuum space. The latter is validated with morphological measurements.


Introduction
Numerical methods are increasingly used for studying the mechanical response and other properties of concrete.The development of computational tools aims at supplementing or even replacing experiments in domains where the latter are difficult to carry out.The response of concrete at very long time scales, typically creep, is but one example.Reproducibility is another main concern, in materials where a large set of aggregate types and sizes are employed.
The use of homogenization tools to tackle full-field micro-mechanical problems is, however, challenging (Dunant et al., 2013).The mechanical and durability properties depend to a large extent on the shape and spatial arrangement of aggregates and pores inside a multi-scale microstructure.Other structure-dependent physical properties include fracture and crack propagation (Stroeven, 2000) or the spatial distribution of highest stressed zones (Escoda et al., 2011b).Thus, a very fine description of the material microstructure is needed for numerical approaches.Numerical methods in general apply to simulated microstructure models or to segmented images of real microstructures, be it concrete (Nagai et al., 1998), mortar (Escoda et al., 2011b) or cement paste (Gallucci et al., 2007;Bullard & Garboczi, 2006;Hain & Wriggers, 2007).
The modeling of cement paste usually combines a cement microstructure and hydration simulation (Thomas et al., 2011;Stark, 2011).Authors have studied the influence of water-to-cement ratio, hydration degree, capillar porosity, particles size distribution or image resolution.In a number of works, spherical shapes (Bentz, 1997;Haecker et al., 2005;Šmilauer & Bažant, 2010;Bentz et al., 1999;Bullard & Garboczi, 2006;Garboczi & Bentz, 1992;Ye et al., 2003;Bishnoi & Scrivener, 2009) are used to model cement (e.g.before hydration).Similar methods are applied to mortar modeling, e.g.Bernard et al. (2008).Other models for cement paste directly address the problem of the cement paste microstructure after hydration.In (Bary et al., 2009), the cement paste is a threephases microstructure made of aggregates and voids in the C-S-H matrix.The authors make use of spherical aggregates or a mix of spherical and prismatic shapes.Like cement aggregates, concrete aggregates are often modeled by discs in 2D (Schlangen & van Mier, 1992;Schutter & Taerwe, 1993;Zaitsev & Wittmann, 1981), or by spheres in 3D (Bazant et al., 1990;Gal & Kryvoruk, 2011;Wriggers & Moftah, 2006).Specifically, the models of Schutter & Taerwe (1993) and of Pedersen et al. (2007) allow a granular skeleton compaction of 1100 kg/m 3 .Other models not based on spheres or discs were also developed.Zaitsev & Wittmann (1981) considered aggregates with polygonal shapes.Similarly, crushed aggregates have been modeled by polygons in 2D (Wang et al., 1999;Kwan et al., 1999).In these models, the number of edges of the polygons as well as their length and orientation are random variables, and aggregates are elongated and scaled before being implanted in decreasing size order.The 2D model of Dequiedt et al. (2001) involves a Voronoï tessellation followed by an erosion and an opening.Voids are modelled in the matrix phase by a Boolean model of spheres whose radius follows an exponential law.Kim & Abu Al-Rub (2011) tested in 2D the influence of the aggregates shape on void nucleation and crack propagation.Discs, hexagons, pentagons, tetragons and arbitrary polygonal shapes have been considered.The influence of the size distribution of aggregates, their volume fraction and the thickness of an interface transition zone has been studied as well.
In 3D, models using ellipsoidal and polyhedra-shaped aggregates have been compared in terms of sphericity and specific surface area (He et al., 2010).Among the shapes the authors tested, octahedra provided for the closest adjustment to X-ray tomography data.Modified ellipsoids with sinusoidal surfaces have been used by Häfner et al. (2006) to model 3D aggregates.A Voronoï tessellation with a BCC distribution of germs have been used by Caballero et al. (2006) to model 3D microstructures.In other works, tomography images of cement and concrete have been used to extract and model aggregates by means of harmonic spherical functions (Garboczi & Bentz, 1992;Garboczi, 2002).Static or dynamic discrete element modeling, which mimics the packing process of particulate materials, are frequently used to model the dispersion of particles in highly-dense materials like concrete (Stroeven et al., 2012a,b).
At the scale considered in this paper, concrete is most often made of a cementitious paste surrounding aggregates and voids.The aim of this study is to provide a model that accounts for two families of morphological criteria concerning aggregates in concrete, namely their size distribution, as given in the formulation, and their distribution in space as described by the covariance of the material estimated on 3D images.In particular the behaviour of the covariance for short distances is related to the ratio of the surface over volume of aggregates, while its long range behaviour accounts for the scale induced by the packing.Our model is based on morphological data obtained from 3D micro-tomography images and from concrete's formulation.Formulation should be interpreted here as sieve curve for concrete technologists.One sample of a concrete material is first introduced in section (2).Its formulation and 3D tomography images are given.The sample's morphological characteristics are derived from the 3D image in section (3).In section (4), a model for concrete microstructure simulation is developed: two types of polyhedra family, Poisson polyhedra and Voronoï polyhedra are considered and compared to the model's aggregates.In section (5) we propose a microstructure model and validate it using morphological measurements.Moreover, an algorithm used to generate Poisson polyhedra simulation is introduced and validated in Appendix (A).This algorithm is "vectorial", i.e. all computations are analytically carried out in the continuous space.The choice of polyhedra is adequate to reproduce the morphology of aggregates obtained by fragmentation (or crushed rock aggregate) as in the present material.Other shapes would be required for other types of materials.We are not considering here details of the shape of polyhedra like statistics about their number of faces or edges, but only their volume distribution, and their overall average shape as contained in the covariance.

Concrete sample
This work is based on a concrete sample that is introduced in this section.Hereafter we present the concrete's formulation granulometry and 3D tomography images.

Morphological data
The aggregates formulation of the concrete sample is divided in three classes of aggregates diameters: sand (0 − 4 mm), fine gravels (4 − 12. 5 mm) and gravels (10 − 20 mm), with relative proportions of 0.41, 0.17 and 0.42, respectively.The three classes of aggregates granulometries are noted G 0−4 , G 4−12.5 and G 10−20 for sand, fine gravels and gravels, respectively.They are given in Fig. (1a).The aggregates full granulometry G(r) is derived according to the relative proportions of each class as: as shown in Fig. (1b).Each of the three functions entering the right-hand term above is given numerically for a finite number of values of r.Linear interpolation has been used to apply the formula as the values of r do not coincide for the three granulometries.We emphasize that the material studied in this work is a simplified model of concrete with a limited number of sieves.The latter is nevertheless useful for simulation purposes.
Hereafter, aggregates with a diameter smaller than 0. 1 mm are incorporated in the cement paste.The latter is modelled as a homogeneous phase.From the mass/volume density of aggregates as well as the volume of the concrete sample, we deduce the volume fraction of the three classes: 25.9%, 12.3% and 30.5%, respectively, for classes 0.1 − 4 mm, 4 − 12. 5 mm and 10 − 20 mm.Hence, the total aggregates volume fraction is 68.7%.

3D Microtomography
A 3D image of a concrete sample is acquired using microtomography technique.This image is statistically representative of the material if the sample is large enough.Concrete is a material with high absorption.To let X-rays cross the sample, a powerful source is needed.This is especially important to obtain contrasted images with a low signal-to-noise ratio.This acquisition was carried out at École des Ponts ParisTech (Marne-la-Vallée, France), on a cylindrical sample of diameter 110 mm.The image is of size 1224 3 voxels, for a resolution of 54 µm/voxel.The 16 bits data image was reconstructed from 1440 projections.A 2D cut of the microtomography image is shown in Fig.
(2).In the next section, we use the tomography image to determine the sample's morphological characteristics.

Morphological characteristics of the concrete sample
In this section, we first correct the tomography image bias regarding non-uniform luminosity and noise.An improved gray-level image is obtained.We then derive some of the concrete's morphological characteristics by measurements carried out on the improved gray-level 3D tomography image.Hereafter we denote by g(x, y, z) the value of the original tomography image along voxel of coordinates x, y and z.The image size is

Luminosity bias correction and noise filtering
As a first improvement of the microtomography image, we correct two types of luminosity bias, characterized by non-uniform luminosity over the sample.For an image of good quality, the vertical variations of illumination should be homogeneous, without any drift as observed on the reconstructed image.Hereafter, we correct the drifts by removing on each horizontal section its average illumination.Without this correction, the aggregates segmentation leads to non-uniform density.The first "vertical" bias is easily noticed when computing the means of the function g along planes perpendicular to the z-axis, as a function of the coordinate z (Fig. 3), i.e.
The z-axis is also the axis along which the sample was rotated during the acquisition process.Luminosity variations of up to 8% are observed.Uncorrected, this drift of the luminosity leads to biased morphological measurements such as anisotropic covariance functions, whereas the material is presumably isotropic.Accurate covariances are in particular necessary for the validation of the morphological modeling of the material.
The second "radial" bias is well-known on microtomography images (Escoda et al., 2011b).The means of g over lines parallel to z is represented as a function of the coordinates x and y in Fig. (4a).Concentric circles appear along the rotation axis.
The vertical and radial luminosity biases outlined above are removed by computing local means and substracting them from the initial image.However, the unwanted influence of well-contrasted voids does not allow one to compute these local means directly from the initial image g, as pointed out in Escoda et al. (2011b).Accordingly, a treatment of the voids is carried out as a preprocessing step.This consists in extracting voids, by a manual segmentation with threshold 18000, and replacing them by the global mean over the whole sample g : The segmentation of the voids used above needs not be an accurate representation of the porous phase.Hereafter, the subscripts 2D , 1D and the superscript ′ are used to indicate the transforms given in equations ( 2), ( 3) and ( 4) resp.The "vertical" and "radial" luminosity biases are then corrected by computing successively A 2D cross-section of the corrected image g 2 is shown in Fig. ( 5).The 2D image (g ′ 1 ) 1D used in the luminosity correction is shown in Fig. (4b).The effect of voids has been reduced compared to image (4a) that represents g 1D .However, effects of the large aggregates are still present.We now remove some of the noise from the corrected g 2 image using a 3D Gaussian filter, following Escoda et al. (2011b).The filter depends on a size parameter R. We consider values of R in the range 1 ≤ R ≤ 6 in voxels.For each value of R, we apply an automatic thresholding on the filtered image.We determine the value of R "visually" as the one giving the "best" automatic thresholding.The value R = 4 is chosen.A 2D cross-section of the filtered image is shown in Fig. (6).Use of a median filter gives similar results, however the Gaussian filter provides smooth boundaries after segmentation, and is therefore more appropriate here.

Morphological measurements
In this section, the granulometry and covariance of the granular phase are measured from the grey-level image.These morphological measurements are used in section (4.2) to derive and validate a random morphological model for concrete.Voids are first "eliminated" so as to obtain measures that are relevant to the aggregates phase.We segment them using an automatic thresholding that maximizes the entropy of two-phases media (Escoda et al., 2011b;Kapur et al., 1985;Sezgin & Sankur, 2004).We then set the porous phase to a constant value, equal to the mean over the cement paste phase.Again, the latter is determined using automatic thresholding (Otsu, 1979;Sezgin & Sankur, 2004).

Granulometry of aggregates
The granulometry is first measured on the grey-level image.The cumulative granulometry G A (r) of a set A (i.e. a binary image) is obtained by: where the set A Fr is the opening (i.e. one erosion followed by a dilation) of A by the structuring element F r .We choose a rhombicuboctaedron of radius r for the structuring element F r , to approximate a ball of radius r.For the grey-level image g(x), the cumulative granulometry is: where g Fr (x) is the opening of g by F r , and E is the expectation of the grey-level image.
The granulometry G(r) measured on the 3D grey-level image, as well as concrete formulation, are given in Fig. (7).The two granulometries are different.First, aggregates of radius larger than 6. 25 mm are not correctly characterized, because the volume of the sample is not large enough to be representative of the material.Second, the image resolution is 0.054 mm, which is not sufficient for characterizing aggregates with radius below this value.This effect is magnified by the low signal-to-noise ratio of the image and the filtering of high frequencies used to suppress noise.As the sensor size of the tomography device is fixed, a compromise should be found between volume size and resolution when acquiring images for multi-scale materials such as concrete.

Covariance of aggregates
The covariance function is a measure of the correlations in the microstructure at varying length scales.It is also a sensitive criterion for identifying concrete models representative of the real microstructure.The centered covariance is defined by: whereas the correlation function is: with D 2 g as the variance over the image of grey-level g(x).The correlation function, shown in Fig. ( 8), is computed on the grey-level image after luminosity bias correction, noise filtering and void removal.This function will be used in section (4.2) for validating the morphological model.
Granulometry measurements from the microtomographic image show that the examined sample is not representative of the material, in terms of size distribution of aggregates.One should generate representative and realistic random models to study the physical behavior of concrete.We develop such a model in the remaining part of this paper.

3D morphological modeling of concrete
In this section, we develop a morphological model for concrete microstructure simulations based on the contracted formulation of the aggregate size distribution.Two types of polyhedra are first considered to model the aggregates shapes: Poisson and Voronoï polyhedra.The experimental aggregate size distribution is then modeled by truncated multi-scale size distribution of each family of polyhedra.A model based on multi-scale Poisson polyhedra is chosen and validated.

Poisson point process
In the homogeneous case, the Poisson point process, of density θ ≥ 0 in the space R 3 is a random point generation process such that the number of points N (K) contained in a compact K follows a Poisson law with expectation N 0 = θµ 3 (K) (with µ 3 denoting the Lebesgue measure, or volume, in R 3 ): For a Poisson point process, the numbers N (K i ) are independent random variables for any family of disjoint compact sets Figure 7: Granulometry measured on the grey-level image after "eliminating" voids (red) vs. concrete formulation (black): cumulative granulometry G(r) (a) and granulometry distributions g(r), as functions of the radius r.

Poisson tessellation
A Poisson tessellation of R 3 is generated by implanting planes.In its homogeneous and isotropic variant, it is a Poisson point process in the space S 3 × R + , where S 3 is the unit sphere and R + is the set of positive real numbers (Matheron, 1971(Matheron, , 1972)).Each plane H of the tessellation is identified by its normal u and distance to the origin r.
The homogeneous isotropic Poisson tessellation model is completely determined by its density λ ≥ 0. The density of planes λ is the Poisson point density induced on any line of the space R 3 .This is measured by computing the density of the number of intersects between the planes and one line.

Parameter assessment for multiscale modeling of aggregates
In this section, models of concrete aggregates are introduced.Both Poisson and Voronoï polyhedra are considered to represent the concrete's granulometry formulation.The formulation of concrete is made of three classes of aggregates: sand, fine gravels and gravels.The characteristic length scales in each class are well separated as shown by the convexity changes in the cumulative granulometry (Fig. 1b).We accordingly model the concrete's formulation granulometry (i.e. its sieve curve) by a combination of three truncated granulometries.

Concrete formulation by multiscale truncated granulometry
Hereafter, the concrete formulation is successively modeled by a combination of truncated granulometries of Poisson polyhedra (section 4.2.2) and of Voronoï polyhedra (section 4.2.3).We denote by G P (r; λ) and G V (r; V m ) the full (non-truncated) cumulative granulometries of Poisson and Voronoï polyhedra, respectively.The former depends on the Poisson planes intensity λ and the latter on the average polyhedra volume V m .Both granulometries are generically noted G M (r; ξ) where M = V or M = P and ξ = λ or ξ = V m is related to a size parameter.The truncated cumulative granulometry G M Λ 1 ,Λ 2 (r; ξ) is the cumulative granulometry of Poisson (or Voronoï) polyhedra with radius comprised between Λ 1 and Λ 2 .We recall that the radius of a polyhedron is the minimum size of the structuring element necessary to remove the object using an opening by reconstruction.Accordingly, the truncated granulometry is given, in terms of the full granulometry, by the conditional law: ).We now model the concrete formulation granulometry, refering to Eq. ( 1) by: where Λ 0 = 0, Λ 3 = 12.5mm consistently with figure (1b).The three granulometries on the right-hand side are truncated in the intervals [Λ i−1 , Λ i ] (1 ≤ i ≤ 3) and the relative weights are as in Eq. (1).We let Λ 1 and Λ 2 vary because the granulometry classes slightly overlap one another.We determine the two truncation parameters Λ 1 and Λ 2 and the three size parameters ξ i (i = 1, 2, 3) by minimizing: where the sum is carried out over all values of r given by the concrete formulation (Fig. 1b).This least-square minimization is carried out numerically using the nonlinear optimization algorithm in Nelder & Mead (1965).We emphasize that the aggregates's sizes in Eq. ( 1) are diameters measured by sieving, whereas Λ 1 and Λ 2 are radii.Hereafter, we use radii instead of diameters when referring to the aggregates's sizes.The latter are easier to work with when using morphological operators.

Model parameters for Poisson polyhedra
The granulometry by opening G P * (r; λ) of a set of Poisson polyhedra of intensity λ is given by (Matheron, 1971(Matheron, , 1972)): We emphasize that the granulometry considered in section (4.2.1) makes use of openings by reconstruction, whereas the above one is valid for a classical openings granulometry.Nevertheless, numerical results carried out in appendix A.2.4 (see also Fig. 25) show that both granulometries are very close to each other so that we assume G P (r; λ) ≈ G P * (r; λ).The expression in ( 13) is found to be minimal (equal to about 0.012) when: The corresponding granulometries are shown in Fig. (10).The values Λ 1 and Λ 2 are estimates of the minimum and maximum radius in the second class, and were computed by the optimization process.They compare well to the diameters given by the formulation (section 2.1), which are respectively 2 and 6.25 mm.

Model parameters for Voronoï polyhedra
Contrary to Poisson polyhedra, no exact theoretical expression is known for the granulometry G V (r; V m ) of Voronoï polyhedra.Nevertheless, Marthinsen (1996) carried out simulations and proposed approximations using a gamma law.We recall the expression of the gamma law of parameters η > 0 and ϑ > 0: where Γ is the Gamma function, i.e.Γ(t) = +∞ 0 dt t x−1 e −t .The expectation of the gamma law is ηϑ.The distribution in number f X of the random variable X = V /V m was fitted by the gamma law (Marthinsen, 1996): where k = 5.56.The granulometry in volume g X corresponding to f X is g X (x) = ̺V f X (x) where ̺ > 0 is chosen so that g X is of mean 1. Making use of dx xf X (x) = 1, we find that ̺ = 1/V m , and: Similarly, the distribution in volume of the random variable V is: For comparison with the concrete formulation, this distribution should be expressed in terms of radius r.We use: where we have assumed V is a sphere.The expression in ( 13) is found to be minimal when: The corresponding granulometries are represented in Fig. (11).The fit is closer to the experimental data with Poisson polyhedra (score 0.012) than with Vornoï polyhedra (score 0.053).This result is consistent with previous models of materials obtained by iterated fragmentation using Poisson polyhedra (Matheron, 1971).Hereafter, we use Poisson polyhedra solely to represent the concrete aggregates.
5 Modeling of the concrete microstructure using Poisson polyhedra We detail hereafter the selection and implantation of Poisson polyhedra following the granulometry fitted in Fig. (10a).Polyhedra are selected from a pre-computed library of Poisson polyhedra.After implantation, the polyhedra do not overlap, as aggregates cannot physically intersect.

Polyhedra selection
We detail hereafter how aggregates are selected from a library of pre-computed Poisson polyhedra.Each polyhedron is stored in vector description in the library together with We first apply an homothethy of factor α to all polyhedra to account for the change of scales between the library and that of the simulation: where L is the size of the simulated microstructure and R i is the polyhedron radius.
The quantity 2π √ 3λ I L represents the average number of Poisson planes intersecting the sphere circumscribed to the volume of the simulation (see appendix A).We retain polyhedra with radius R ′ i in the range Λ 1 < R ′ i < Λ 2 .To achieve the volume fraction f I , we choose randomly in the retained Poisson polyhedra n p aggregates with: where K T 0 is the average Poisson polyhedra volume in the truncated granulometry.In turn, we compute K T 0 numerically using the formula of Miles-Lantuéjoul, to correct for border effects (see appendix A): where α 3 accounts for the scale change between the library and the scale of the simulation, B i and V i are the bounding box and the volume of the polyhedron i, respectively.
D is the simulated volume and I is the set of polyhedra in the truncated granulometry.The formula for the probabilty P [B i ⊂ D] is given in appendix (A).

Algorithm for polyhedra implantation
Polyhedra implantation is carried out as detailed below.No overlap occurs but contacts are allowed.We suppose that the set of Poisson polyhedra to be implanted is given and follows the volume fraction of each compacted class of aggregates in concrete, i.e. 25.9%, 12.3% and 30.5% of polyhedra for classes 0.05 − 2. 4 mm, 2.4 − 6. 15 mm and 6.15 − 12. 5 mm, respectively (see section 2.1).We first order polyhedra by decreasing volume, irrespective of the classes they belong to.Large aggregates are implanted first to maximize the total aggregates volume fraction that the method achieves.Use of the volume of polyhedra as the sorting criterion is appropriate in the context of convex polyhedra with moderate aspect ratios.Prior to the implantation of a polyhedron, we generate two vectors: a position x uniformly distributed in the 3D space and a direction v uniformly distributed in the unit sphere.The polyhedron is initially implanted in x or equivalently, translated by a vector x from the origin.The packing is carried out using an attraction-repulsion mechanism through a translation of the polyhedron in the direction v. Intersection between the current polyhedron and others already-implanted polyhedra are tested.At this point, two possibilities arise.In the case of an intersection, a repulsion mechanism is used.The polyhedron is iteratively translated in the direction v by step of size ̟, until there is no intersection anymore.When one-quarter of the size of the image has been scanned, a new position x and a new direction v are generated and used for this polyhedron.In the case of no intersection, an attraction mechanism is used.The polyhedron is translated in the direction v by steps ̟ until there is a non empty intersection.When there is intersection, the polyhedron is moved once in the opposite direction, with a translation length ̟.The repulsion mechanism ensures there is no overlap between implanted polyhedra.The attraction mechanism provides higher maximum volume fractions of polyhedra in the microstructure.To avoid creating clusters containing large polyhedra, which are implanted first, the first ten polyhedra are implanted without attraction: only the repulsion mechanism is used.
The above procedure does not guarantee that a polyhedron will be implanted.When a large number N max of initial positions and directions followed by the repulsion mechanism have been generated for a given polyhedron, we consider that the polyhedron cannot be implanted and that the desired volume fraction is too high to be attained.We choose N max = 1000.
Translations and scalings are performed in the same manner as for the implantation of polyhedra in a Boolean model, detailed in appendix (A.2.5).Although the polyhedra are defined vectorially in the continuous space, testing for polyhedra intersections vectorially would be too cumbersome.In the procedure outlined above, we test intersections of Poisson polyhedra using discretizations on a grid of voxels.This binarization is a convenient intermediate step in the computation.In the end, the method nevertheless allows one to simulate vectorial microstructures of polyhedra packings.Finally, we also choose to generate periodic microstructures for later use with Fourier-based methods to compute concrete's physical properties (Willot et al., 2014).

Maximum filling ratio
Hereafter we give the maximum filling ratio, or aggregates volume fraction, our method achieves, for the truncated and non-truncated Poisson polyhedra granulometries.For a given granulometry, the target volume fraction is increased until the method does not provide a valid microstructure, i.e. one polyhedron could not be implanted.We fix the size of the microstructure to 500 3 voxels, the step size to ̟ = 1 voxel.We also fix the resolution so that the number of implanted polyhedra should be in the order of 1000, as a compromise between representativeness and computational efficiency.Equation ( 22) provides the resolution in the microstructure in mm per voxel.The filling ratio is augmented by increasing the number of aggregates.For each target volume fraction, 10 microstructures are generated.A filling ratio is attainable if at least one of the 10 images was successfully generated at the given volume fraction.
Results for the compaction tests and the percentage of polyhedra retained when truncating the granulometry are given in Tab. ( 1).The highest aggregates density is attained by the non-truncated granulometry, whereas, as expected, narrow granulometries provides the lowest maximal density.These figures give a first indication of the volume fractions of aggregates of the three classes that we will be able to implant in the microstructure.Suppose that polyhedra from the two truncated granulometries with the largest radius have been implanted with the volume fractions corresponding to the concrete formulation.According to data in section 2.1, a volume space of 57.2% is left in the material for the cement paste and for the third truncated Poisson polyhedra granulometry.The concrete's formulation requires a volume fraction of 25.9% for the latter, so that 45% of the space left has to be filled with the smallest polyhedra.This is slightly larger than the maximum volume fraction of 42% which can be attained in an empty cube.We therefore conclude that it will not be possible to fully respect the concrete's formulation as far as the smallest aggregates are concerned.Nevertheless the difference is small.This is investigated in the following section.

Microstructure generation
We generate a microstructure containing 1600 3 voxels according to the parameters given in Tab.
(2).The number of aggregates is about 10 6 .The volume fractions of aggregates from the two truncated granulometries at the largest scales are that of the concrete's formulation.Smaller aggregates from the third truncated granulometry are implanted with a volume fraction of 24%.As expected, this is slightly less than concrete's formulation.Due to the large number of contact points between implanted polyhedra, the aggregates percolate in the simulated microstructure.Percolation has strong effects on the overall elastic properties of the composite (Willot & Jeulin, 2009).When the aggregates are much stiffer than cement paste, the effective elastic moduli are overestimated.We accordingly disconnect aggregates by disconnecting by applying a 3D watershed algorithm (Beucher & Meyer, 1992) to the inverse of the distance function inside the polyhedra phase (Escoda et al., 2011b).Polyhedra labels are used as markers.The impact on the volume fraction is small: the aggregates volume fraction is down from f t ≈ 65.2% to f LP E t ≈ 64.3% after disconnection in the three-scales microstructure.This is in contrast with materials generated on smaller grids.For a microstructure of size 512 3 voxels with the same Poisson intensity as in Tab. ( 2), but resolution 0. 3 mm/voxel, the aggregates volume fraction decreases from 0.657 to 0.532 after applying the watershed algorithm.Not surprisingly, fine discretization is necessary for generating non-percolating aggregates following concrete's formulation.

Model validation using covariance and granulometries measurements
We first check the influence of the step parameter ̟ on the implantation algorithm with regards to covariances.We consider microstructures containing two classes of polyhedra with the largest sizes and varying step sizes ̟ = 1, 2, 4 voxels (Tab.3).A 2D section of the microstructure with ̟ = 2 is represented in  (3) with step size ̟ = 2 .
We now compare the centered reduced covariance C c,r (h) of the three-scales model in Tab. ( 2) with the tomography image (Fig. 15).The characteristic lengths of these two covariance graphs are identical, however the curves are very different in the range h < 10 mm.This is presumably an effect of the low resolution of the tomography image and of the absence of small aggregates, in the class 2 − 6.25 mm, as explained in section (3.2.1).These small aggregates are taken into account in the simulated microstructure making use of the formulation.
As a criterion, the covariance function of the tomography image is not relevant given the very large differences between the smallest and largest aggregates.We therefore did not investigate the covariance properties of a model based on Voronoï polyhedra.
Finally, we now investigate the granulometry of our generated model and compare it to that of the tomography image (Fig. 16).Opening by reconstruction granulometry is measured on the three-scales microstructure.In terms of the cumulative granulometry, the largest difference between the simulated granulometry and that of the theoretical one is found to be about 8% at r = 6.5 mm.We interpret the latter as size effects at the highest length scale, which are due to the limited representativity of the simulated multi-scale microstructures.In that respect, the difference should be compared with fluctuations in the granulometry of the generated model, which are about 4% when r is large.

Conclusion
The procedure presented in this work allows one to simulate concrete microstructures using experimental granulometry data.Our model takes into account the granulometry of concrete aggregates in the range 0.1 − 20 mm, reproducing the size distribution of a simplified model material.The concrete's formulation granulometry is approached by  combinations of three families of polyhedra, representing sand, fine gravels and coarse gravels.Poisson and Voronoï polyhedra were emphasized for this purpose, with the best results achieved with the first type.
We developped and validated an algorithm for generating vectorially and implanting Poisson polyhedra in a fine-grained mortar matrix.The model parameters are derived from the granulometry and volume fraction of aggregates specified by the concrete formulation.They are efficiently determined by optimizating on analytical formula.More generally, our model could be applied to other granular media of similarly-shaped grains.The simulated microstructures allows one to predict concrete's local and effective properties (Escoda et al., 2011b).In that respect, comparisons between our Poisson-polyhedra model and binarized microtomography images should be carried out in terms not only of the effective but also on the local elastic response Escoda (2012).A Appendix: Implementation and validation of an algorithm to simulate Poisson polyhedra We give in this appendix a vector implementation of an algorithm that generates Poisson polyhedra.This extends the presentation given in Escoda et al. (2011a).The vector implementation outlined here consists in an analytical definition of the polyhedra contrary to bitmap implementations in which images are defined on a grid of voxels.
Vector algorithms allow for faster generation and have lower memory cost, both in terms of RAM when generating the polyhedra and in terms of hard drive memory for storing.Last, a polyhedron defined vectorially may be binarized at arbitrary resolutions.
Our algorithm is validated by comparing theoretical expressions of various morphological measurements (Matheron, 1975) and measurements carried out on simulations.This algorithm is used to create a library of vector Poisson polyhedra.An open-source version of the algorithm and code will be made available as part of EDF's Material Ageing Platform (http://materialsageingplatform.org).A vectorial library of Poisson polyhedra will be shared at http://cmm.ensmp.fr/∼willot.

A.1.1 Space tessellation
A Poisson tessellation of a cubic subdomain D of size L included in the space R 3 is first generated, following section (4.1).We limit ourselves to a Poisson points process in S 3 × [0; r 1 ], where r 1 = √ 3L/2 is the radius of the circumscribed sphere to the cubic domain D, to avoid generating planes which do not intersect D. The expectation N of the number of Poisson planes intersecting the circumscribed sphere is proportional to λ and to the measure of the set S 3 × [0; r 1 ]: We implant N P planes, where N P follows a Poisson distribution (Eq.10) with expectation N above.The planes (P i ) i≤N P are defined by their normal vectors u i , taken uniformly on the sphere S 3 and by the distance r i that is uniformly distributed over [0, √ 3L/2].The uniform distribution on the unit sphere S 3 is obtained by the method described in Knuth (1969): we generate a random vector (x 1 , x 2 , x 3 ), where the x i are random variables following a standard normal distribution; the normalized variable x/ x is uniformly distributed on the unit sphere S 3 .

A.1.2 Polyhedra labeling
To select polyhedra uniformly random in number, by contrast to a selection in volume, we label all polyhedra in the tessellation.Additionaly, only polyhedra that are not cut by the boundary of D are labeled.The tessellation is defined as a set of planes (Π i ) i≤Np of Cartesian equations a i x + b i y + c i z + d i = 0.Each polyhedron is completely defined by its location in one of two half-spaces delimited by each plane.Thus a polyhedron is labelled as a list of signs (s i ) i≤N P such that all points (x, y, z) in the polyhedron verify: We now consider all straight lines (∆ i,j ) i,j≤N P at the intersection of pairs of planes Π i and Π j , as illustrated in figure (17) in 2D.We first eliminate lines which do not intersect the domain D of the image.Consider now a particular line ∆ i,j among the remaing ones and its intersections with all planes in the tessellation (Π kc ) kc≤Np .We select the M − 1 intersecting points (I k ) 1≤k≤M −1 that lie inside D, ordered by their x coordinates so that x k < x k+1 .By convention, we set I 0 and I M as the intersections of ∆ i,j with D, so that ∆ i,j is cut into M segments.Each one of these segments is the edge of four polyhedra in the tessellation.The values of (s k ) k =i,k =j are common to all four polyhedra.They are equal to the sign of a i x + b i y + c i z + d i at any point in the segment [I K ; I K+1 ], for instance (I k + I k+1 )/2.The signs (s i , s j ) take on the four possible values (±1, ±1) for each of the four polyhedra.
Carefulness is required to determine incomplete Poisson polyhedra, i.e.Poisson polyhedra that intersects D. We also wish to determine the list of planes that form the faces of each polyhedron and its bounding box, useful for binarization.To this aim, we maintain a flag indicating wether each polyhedron is a "complete" or "incomplete" one.This flag is updated during the iterations process as follows.Each polyhedron is added if not already present.When considering the extremal segments [I 0 ; I 1 ] and [I M −1 ; I M ], the flag is set to "incomplete" irrespective of its previous value.When considering interior segments [I k ; I k+1 ] (0 < k < M − 1), the flag is unchanged if the polyhedron is already present, or set to "complete" if the polyhedron is new.When a polyheron with a "complete" flag is updated or created, we additionally update its bounding box and set of faces.The bounding box of a polyhedron is the smallest rectangular parallelepiped circumscribed to the polyhedron.To simplify, only bounding boxes with faces parallel to the reference axes are considered so that the faces equations are of the form x = x min and x = x max (likewise for y and z).When iterating over a given segment [I k ; I k+1 ] we update the values of x min by x min ← min(x min , x I k , x I k+1 ) and similarly for and x max and the other coordinates.Finally, the set of faces are updated by adding (if not already present) the planes Π i and Π j as well as the two planes that correspond to the intersection points I K and I K+1 . is the probability that the bounding box B of an object is fully enclosed in D (Serra, 1982).It is given by: where L is the size of the cube D and L x,y,z are the dimensions of B. This is schematically illustrated in Fig. (18).

A.1.4 Binarization of a polyhedron
We use the list of signs to binarize a given Poisson polyhedron on a 3D image grid.To speed up computations, only signs related to planes that define the polyhedron faces are checked.Furthermore, we make use of the convexity property of Poisson polyhedra, as show in Fig.  A.2 Validation of the algorithm and polyhedra library

A.2.1 Geometrical covariogram
We recall the expression for the geometrical covariogram K(h) of a set A: where µ 3 is the Lebesgue measure on R 3 , and A −h is the translation of A by the vector h, as shown in Fig. (21).For Poisson polyhedra of Poisson planes intensity λ: where K(0) is the volume expectation and h is the norm of h.We now compute numerically the geometrical covariogram on polyhedra that were produced by the method described in the previous section.Various values of the expected number of Poisson planes N = 80, 100, 150, 200, 300 are investigated, and the geometrical covariograms are computed for vectors h in the x direction, after discretization of the domain D on a grid of size 256 3 voxels.The geometric covariogram is computed for N t = 10000 realizations (resp.N t = 5000) when N = 80 or 100 (resp.N = 150, 200, 300).For each value of N , a family of polyhedra (A i ) i≤Nt is thus generated.The geometrical covariogram is estimated by averaging: We also compute:      Jeulin (1991)).We note A the resulting set.An example of a 2D cross-section of A is shown in Fig. (23).Its covariance reads: where r(h) = K(h)/K(0) is the normalized geometrical covariogram given in Eqs. ( 27) and ( 28).The volume fraction of the set A is equal to f .The covariance in ( 31) is estimated in the x-direction as: with K D (h) = L 2 (L−h).As many as N A = 5×10 4 random points uniformly distributed in a domain D are used to estimate K A∩D (h) in vector images, without binarization.It is estimated by: where h is a vector of norm h oriented along the direction x, and g(x) = 1 (resp.

A.2.3 Poisson polyhedra library
We generate a library of 5 × 10 5 polyhedra using a Poisson intensity N = 200.The equations of each polyhedron face and its bounding box are included for easier binarization, as well as its volume and weight for Miles-Lantuéjoul correction.Implantation of Poisson polyhedra with arbitrary intensity λ requires the rescaling given in section (5.1).This is computed analytically.For instance the plane Π i defined by a i x+b i y+c i z+d i = 0 is implanted in point (x c , y c , z c ) and rescaled by α by changing d i into: (34) The labels (s i ) remain unchanged after geometric transformations on the planes.

A.2.4 Granulometry of polyhedra of the library
We now validate the opening granulometry of Poisson polyhedra generated in the library.The radius of an object is the minimum value of R necessary to remove the object by opening with a sphere of radius R. Equivalently, the radius of the object A in the opening granulometry is the radius of the largest sphere inscribed therein.The radius R and center x c of such a sphere are called the Chebyshev radius and center.For a given polyhedron, they are analytically determined by the convex optimization problem: where V i is a vector of dimension 3 such as V i = (−s i a i , −s i b i , −s i c i ) and with v i given by v i = s i d i .We solve the problem above numerically using the convex programming software CVX (Grant & Boyd, 2011, 2008).We then compute the radii histogram over polyedra in the library.We weight the histogram using the volume, to obtain a granulometry in volume, instead of a granulometry in number, and by the Miles-Lantuéjoul correction.This granulometry is compared to the theoretical Poisson polyhedra granulometry given by (Matheron, 1972): We emphasize that the expression above is the theoretical expression for a (classical) opening granulometry, not an opening by reconstruction as computed by the Chebyshev measure.Nevertheless, a comparison between the two granulometries, represented in Fig. ( 25), show that the difference between the two granulometries is negligible.We add the Chebyshev radius of each polyhedron to the polyhedra library.

A.2.5 Covariance of the Boolean model of Poisson polyhedra
In this appendix, we validate the algorithm for Poisson polyhedra generation by measuring the covariance of Poisson polyhedra Boolean media.In the Boolean model considered here, polyhedra are implanted on primary grains following a Poisson point process of density θ.Polyhedra are allowed to interpenetrate, as represented in Fig. (26).The volume fraction f and the covariance Q(h) of the complementary set are given by, respectively (Jeulin, 1991;Quenec'H et al., 1992): where K 0 is the mean volume and r(h) is the normalized geometrical covariogram of the Poisson polyhedra.We estimate the mean volume V 0 of Polyhedra in the library, after applying the Miles-Lantuéjoul correction: where B j is the bounding box and V j is the volume of the polyhedron j in the library.The average volume of implanted polyhedra is then estimated as K 0 = α 3 V 0 , where α is the rescaling factor.We do not use hereafter the theoretical average volume for Poisson polyhedra (Eq.28), because of the bias in the average volumes of generated polyhedra, as explained in section (A.2.1).This bias is a consequence of the selection of random intact polyhedra in a finite domain.We measure the covariance of the complementary set of a Poisson polyhedra Boolean media with intensity λ = 0.045 and volume fraction f = 0.2.The 3D image is discretized on 500 3 voxels.The measured covariance and theoretical expression are compared in Fig. ( 27).Excellent agreement is found between the two.For values of λ smaller than λ = 0.045, the 3D image is not representative and the measured covariance does not follow the theoretical expression (not shown).

A.3 Conclusion
A library of Poisson polyhedra useful for generating random media has been introduced.This library includes the planes defining the faces and edges of the polyhedra, the polyhedra bounding boxes, weights useful for the Miles-Lantuéjoul correction, the volume and Chebyshev radius.We validated the algorithm for the vector simulation of Poisson polyhedra described in section (A.1) in section (A.2) using various morphological measurements.The vector generation of polyhedra allows to obtain an algorithm with low memory requirements.The CPU time T needed for running the algorithm scales as N 3 log(N ) where N is the expected number of planes in the tessellation inside of the domain D of generation.This scaling law is close to the average number of Poisson polyhedra in the tessellation (equations 24 and 28): The generation of a Poisson polyhedron typically requires 1 min 30 s when N = 200.on a computer using Intel Xeon X5650 (2.66 GHz) processors with one core.

Figure 1 :
Figure 1: Granulometries G 0−4 (r), G 4−12.5 (r) and G 10−20 (r) of the three classes of aggregates in the material formulation as a function of the radius r (a) and full granulometry (b) obtained by linear combination of granulometries given in (a).

Figure 3 :
Figure 3: "Vertical" luminosity bias in the z direction, shown by the vertical profiles g 2D (z).

Figure 5 :
Figure 5: 2D cross-section of concrete tomography image g 2 , after correction of the two luminosity drifts.

Figure 6 :
Figure 6: 2D cross-section of the concrete tomography image after luminosity correction and Gaussian filtering of size R = 4.

Figure 8 :Figure 9 :
Figure 8: Centered reduced covariances C c,r (h) (or correlation function) of the aggregates measured on the concrete 3D grey-level image.

Figure 10 :
Figure 10: Concrete formulation (black) and granulometry of the optimized model of Poisson polyhedra (red) based on three truncated granulometries (a).The model parameters, given in (15), have been optimized.The three truncated granulometries are shown separately in (b).

Figure 11 :
Figure 11: Concrete formulation (black) and granulometry of the optimized model of Voronoï polyhedra (red) based on three truncated granulometries (a).The model parameters, given in (20), have been optimized.The three truncated granulometries are shown separately in (b).

Figure 12 :
Figure 12: 2D cross-section of a 3D microstructure of multiscale Poisson polyhedra (in grey) embedded in a cement paste (red), generated according to the parameters given in Tab. 2 (a).Same microstructure after disconnection with the watershed algorithm with cement paste in black and polyhedra in white (b).
Fig. (13), and the covariances of each model are shown in Fig. (14).The three covariances are very close to each other.Additionaly, we verify that the aggregates percolate in the three microstructures, irrespective of the value of ̟.Accordingly, the step size parameter does not influence the covariance, and can not be used to control percolation.Class (mm) 6.15 − 12.5 2.4 − 6.15 λ I (mm −1 ) Parameters used for generating a two-scales microstructure: Poisson intensity λ I (row 2), aggregates volume fractions f I in the two classes (row 3), resolution, image size and step sizes (rows 4-6 resp.).

Figure 13 :
Figure 13: 2D cross-section of the microstructure corresponding to the parameters given in Tab.(3) with step size ̟ = 2 .

Figure 14 :
Figure 14: Centered reduced covariance C c,r (h) of the three microstructures detailed in Tab.(3).

Figure 15 :
Figure 15: Centered reduced covariance C c,r (h) of the three-scales microstructure given in Tab.(2) (black) and of the micro-tomography concrete image (red).

Figure 17 :
Figure 17: Labeling of polygons generated in a 2D Poisson tessellation.The algorithm iterates over the Poisson straight lines.When iterating on the straight line ∆, ∆ is cut into segments [I K ; I K+1 ] corresponding to edges of polygons in the tessellation.Among these polygons, polygons that are represented in green are complete polygons, while polygons that are represented in blue are cut polygons.Hashed green polygons are cut by the boundary of the domain and will be, or have been placed in the list of cut polygons during iterations on dotted straight lines.
(19).Along a given segment in the bounding box, we determine the two points at the intersection of the line and polyhedron.The interior of the polyhedron are the points in-between these two points.Examples of some of the library polyhedra are represented in Fig.(20).

Figure 18 :
Figure 18: Miles-Lantueéjoul correction: probability for an object to have its bounding box B entirely included in a cubic domain D.

Figure 19 :
Figure19: Binarization of a polyhedron on a 3D grid, making use of the convexity of the Poisson polyhedra.The polyhedron is described by a set of lines in a given direction.Along each line, the interior points in the polyhedron (red) are determined using the points at the intersection of the polyhedron faces and the line (red).

Figure 21 :
Figure 21: Measure of the geometrical covariogram K(h) as the volume of A ∩ A −h , shown in gray for a given vector h.

Figure 22 :
Figure 22: Theoretical value and numerical estimates of the geometrical covariogram K(h) and K corr (h) versus h, in log-lin scale.In the numerical simulations, the number of Poisson planes is N = 200 and 5, 000 realizations were averaged on.
g(x) = 0) if x ∈ A (resp.x ∈ A).Additionaly, we average the covariance on 12 configurations of 3D Poisson mosaic with N = 200 and f = 0.25.Results, shown in Fig. (24), validate the generation of the tessellation of Poisson planes in our algorithm.

Figure 23 :
Figure 23: 2D cross-section of a 3D binary Poisson mosaic of size 256 3 and parameters N = 200 and f = 0.25.

Figure 24 :
Figure 24: Comparison between the theoretical and numerical estimate of the covariance C(h) of a 3D binary Poisson mosaic of parameters N = 200 and f = 0.25.

Figure 26 :
Figure 26: 2D cross-section of a Poisson polyhedra Boolean media with intensity λ = 0.045 and volume fraction f = 0.2.The image is discretized in a 500 3 voxels grid.

Figure 27 :
Figure 27: Covariance Q(h) of the complementary set of a Poisson polyhedra Boolean model with intensity λ = 0.045 and volume fraction f = 0.2: comparison between measurement and theoretical expression.

Table 1 :
Maximum filling ratio for the three truncated Poisson polyhedra granulometries and for the full (non-truncated) Poisson polyhedra granulometry.The percentage of polyhedra retained after truncation, by comparison with the full granulometry, is given in row 3.

Table 2 :
The resulting overall aggregates volume fraction is f t ≈ 65.2%, not far from the formulation (68.7%).2Dcross-sections of the microstructure are represented in Fig.(12).
Parameters used for generating a three-scales concrete microstructure: Poisson intensity λ I (row 2), aggregates volume fractions f I in the three classes (row 3), resolution, image size and step size (rows 4-6 resp.).The resulting aggregates volume fraction is f t ≈ 65.2% and f LP E t ≈ 64.3% after polyhedra disconnection.

Table 4 :
Relative errors on the measure K corr (0) of the average volume of Poisson polyhedra with Miles-Lantuéjoul correction (row 2) for increasing values of the average number of planes N (row 1) used in the simulations of the Poisson tessellation.The average number L/L x is given on row 3.A.2.2 Covariance of a Poisson mosaicWe now validate the Poisson tessellation using a covariance function.A binary Poisson mosaic is generated from a Poisson tessellation by selecting randomly each Poisson polyhedron with probability f (