# 3D morphological modeling of concrete using multiscale Poisson polyhedra

Julie Escoda, Dominique Jeulin, François Willot, Charles Toulemonde

► To cite this version:

Julie Escoda, Dominique Jeulin, François Willot, Charles Toulemonde. 3D morphological modeling of concrete using multiscale Poisson polyhedra. Journal of Microscopy, 2015, 258 (1), pp.31-48. 10.1111/jmi.12213. hal-01096766

HAL Id: hal-01096766

https://minesparis-psl.hal.science/hal-01096766

Submitted on 18 Dec 2014

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# 3D morphological modeling of concrete using multiscale Poisson polyhedra

J. Escoda(1)*, D. Jeulin(1), F. Willot(1), C. Toulemonde(2)

(1) Centre de Morphologie Mathématique, Mathématiques et Systèmes, Mines ParisTech, Fontainebleau, France

(2) Département Matériaux et Mécanique des Composants, Électricité de France, Moret-sur-Loing, France

*Corresponding author.

## Abstract

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 different scales. An algorithm is introduced to generate polyhedra packings in the continuum space. The latter is validated with morphological measurements.

**Keywords** Random microstructure models; concrete; multiscale materials; granulometry

## 1 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 mate-
rial microstructure is needed for numerical approaches. Numerical methods in general
apply to simulated microstructure models or to segmented images of real microstruc-
tures, 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 hydra-
tion 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 three-
phases 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 ce-
ment 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 gran-
ular skeleton compaction of 1100 kg/m³. 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, pen-
tagons, 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 com-
pared 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 to-
mography 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.

# 2 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.

## 2.1 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:

$$G(r) = 0.41G_{0-4}(r) + 0.17G_{4-12.5}(r) + 0.42G_{10-20}(r), \quad (1)$$

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%.

<figure><img src="image_2.png" /><figcaption>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).</figcaption></figure>

## 2.2 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³ 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.

# 3 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
<figure><img src="image_3.png" /><figcaption>Figure 2: 2D cross-section of a 3D tomography image of concrete, containing 1224³ voxels. The image resolution is 54 µm/voxel.</figcaption></figure>

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 $L_x \times L_y \times L_z$ in voxels.

## 3.1 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.

$$g_{2D}(z) = \langle g(x, y, z) \rangle_{x,y} = \frac{1}{L_x L_y} \int dx dy g(x, y, z) \quad (2)$$

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$

$$g_{1D}(x, y) = \langle g(x, y, z) \rangle_z = \frac{1}{L_z} \int dz g(x, y, z) \quad (3)$$

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 subtracting 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⟩:

$$g'(x, y, z) = \begin{cases} g(x, y, z) & \text{if } g(x, y, z) > 18000 \\ \langle g \rangle & \text{otherwise.} \end{cases} \qquad (4)$$

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

$$g_1(x, y, z) = g(x, y, z) - (g')_{2D}(z), \quad g_2(x, y, z) = g_1(x, y, z) - (g'_1)_{1D}(z) \quad (5)$$

A 2D cross-section of the corrected image g2 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 g1D. However, effects of the large aggregates are still present.

<figure><img src="image_4.png" /><figcaption>Figure 3: “Vertical” luminosity bias in the z direction, shown by the vertical profiles g2D(z).</figcaption></figure>

We now remove some of the noise from the corrected g2 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.
<figure><img src="image_5.png" /><figcaption>Figure 4: “Radial” luminosity bias, shown by the means $g_{1D}$ (a) and, after void correction, $(g'_{1})_{1D}$ (b).</figcaption></figure>

## 3.2 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).

### 3.2.1 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:

$$G_A(r) = \frac{P\{\mathbf{x} \in A\} - P\{\mathbf{x} \in A_{F_r}\}}{P\{\mathbf{x} \in A\}}, \qquad (6)$$

where the set $A_{F_r}$ 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(\mathbf{x})$, the cumulative granulometry is:

$$G(r) = \frac{\mathbb{E}[g(\mathbf{x})] - \mathbb{E}[g_{F_r}(\mathbf{x})]}{\mathbb{E}[g(\mathbf{x})]}, \qquad (7)$$
<figure><img src="image_6.png" /><figcaption>Figure 5: 2D cross-section of concrete tomography image $g_2$, after correction of the two luminosity drifts.</figcaption></figure>

where $g_{F_r}(\mathbf{x})$ is the opening of $g$ by $F_r$, and $\mathbb{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.

### 3.2.2 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:

$$C_c(\mathbf{h}) = \mathbb{E}[g(\mathbf{x})g(\mathbf{x} + \mathbf{h})] - \mathbb{E}^2[g(\mathbf{x})], \quad (8)$$

whereas the correlation function is:

$$C_{c,r}(\mathbf{h}) = \frac{C_c(\mathbf{h})}{D_g^2} = \frac{C_c(\mathbf{h})}{C_c(\mathbf{0})} \quad (9)$$

with $D_g^2$ as the variance over the image of grey-level $g(\mathbf{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
<figure><img src="image_7.png" /><figcaption>Figure 6: 2D cross-section of the concrete tomography image after luminosity correction and Gaussian filtering of size R = 4.</figcaption></figure>

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.

# 4 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.

## 4.1 Definitions

### 4.1.1 Poisson point process

In the homogeneous case, the Poisson point process, of density $\theta \ge 0$ in the space $\mathbb{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 = \theta\mu_3(K)$ (with $\mu_3$ denoting the Lebesgue measure, or volume, in $\mathbb{R}^3$):

$$P\{N(K) = k\} = \frac{e^{-N_0} N_0^k}{k!} \qquad (10)$$

For a Poisson point process, the numbers $N(K_i)$ are independent random variables for any family of disjoint compact sets $(K_i)_{i \in I}$, $I \subset \mathbb{N}$.
<figure><img src="image_8.png" /><figcaption>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$.</figcaption></figure>

### 4.1.2 Poisson tessellation

A Poisson tessellation of $\mathbb{R}^3$ is generated by implanting planes. In its homogeneous and isotropic variant, it is a Poisson point process in the space $S_3 \times \mathbb{R}^+$, where $S_3$ is the unit sphere and $\mathbb{R}^+$ is the set of positive real numbers (Matheron, 1971, 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 $\lambda \ge 0$. The density of planes $\lambda$ is the Poisson point density induced on any line of the space $\mathbb{R}^3$. This is measured by computing the density of the number of intersects between the planes and one line.

<figure><img src="image_9.png" /><figcaption>Figure 8: Centered reduced covariances $C_{c,r}(h)$ (or correlation function) of the aggregates measured on the concrete 3D grey-level image.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Figure 9: Voronoï tessellation (in red) with underlying Poisson point process shown as gray spheres (a); Voronoï polyhedron (b).</figcaption></figure>

### 4.1.3 Poisson vs. Voronoï polyhedra

Poisson polyhedra are the complementary set of the Poisson tessellation. A Poisson polyhedron is generated by randomly selecting one polyhedron enclosed by Poisson planes. Polyhedra are selected uniformly random in number, i.e. with the same probability. This is in contrast with a selection in volume where the probability assigned to each polyhedron is proportional to its volume. The latter is obtained by choosing the polyhedron containing the origin. Contrary to the selection in volume, a selection in number requires one to label all polyhedra present in the tessellation. Note that a new tessellation is generated for generating each polyhedron, so that polyhedra are independent.

A Voronoï tessellation of the space $\mathbb{R}^3$ is obtained from a Poisson point process in $\mathbb{R}^3$ (Gilbert, 1962). The space is divided into cells corresponding to the Poisson points influence zones. More exactly, each cell consists of all points closer to one Poisson point than to any other points. The resulting tessellation, shown in Fig. (9) consists of convex polyhedra and each cell is obtained by the intersection of half-spaces.

## 4.2 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.

### 4.2.1 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; \lambda)$ 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 $\lambda$ and the latter on the average polyhedra volume $V_m$. Both granulometries are generically noted $G^M(r; \xi)$ where $M = V$ or $M = P$ and $\xi = \lambda$ or $\xi = V_m$ is related to a size parameter. The truncated cumulative granulometry $G^M_{\Lambda_1, \Lambda_2}(r; \xi)$ is the cumulative granulometry of Poisson (or Voronoï) polyhedra with radius comprised between $\Lambda_1$ and $\Lambda_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:

$$G_{\Lambda_1, \Lambda_2}(r; \xi) = \frac{G(r; \xi) - G(\Lambda_1; \xi)}{G(\Lambda_2; \xi) - G(\Lambda_1; \xi)} \quad \text{if} \quad \Lambda_1 \le r \le \Lambda_2 \quad (11)$$

and $G_{\Lambda_1, \Lambda_2}(r; \xi) \equiv 0$ (resp. 1) if $r \le \Lambda_1$ (resp. $r \ge \Lambda_2$). We now model the concrete formulation granulometry, refering to Eq. (1) by:

$$G_{\Lambda_1, \Lambda_2}^M(r; \xi_1, \xi_2, \xi_3) = 0.41 G_{\Lambda_0, \Lambda_1}^M(r; \xi_1) + 0.17 G_{\Lambda_1, \Lambda_2}^M(r; \xi_2) + 0.42 G_{\Lambda_2, \Lambda_3}^M(r; \xi_3), \quad (12)$$

where $\Lambda_0 = 0$, $\Lambda_3 = 12.5\text{mm}$ consistently with figure (1b). The three granulometries on the right-hand side are truncated in the intervals $[\Lambda_{i-1}, \Lambda_i]$ ($1 \le i \le 3$) and the relative weights are as in Eq. (1). We let $\Lambda_1$ and $\Lambda_2$ vary because the granulometry classes slightly overlap one another. We determine the two truncation parameters $\Lambda_1$ and $\Lambda_2$ and the three size parameters $\xi_i$ ($i = 1, 2, 3$) by minimizing:

$$\sum_r [G_{\Lambda_1, \Lambda_2}^M(r; \xi_1, \xi_2, \xi_3) - G(r)]^2, \quad (13)$$

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 $\Lambda_1$ and $\Lambda_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.

### 4.2.2 Model parameters for Poisson polyhedra

The granulometry by opening $G_*^P(r; \lambda)$ of a set of Poisson polyhedra of intensity $\lambda$ is given by (Matheron, 1971, 1972):

$$G_*^P(r; \lambda) = 1 - \left( 1 + 4\pi\lambda r + \frac{\pi^4}{6}\lambda^2 r^2 + \frac{2}{9}\pi^5\lambda^3 r^3 \right) e^{-4\pi\lambda r}. \quad (14)$$
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; \lambda) \approx G_*^P(r; \lambda)$. The expression in (13) is found to be minimal (equal to about 0.012) when:

$$ \begin{aligned} \lambda_1 &= 0.518 \, \text{mm}^{-1}, \quad \lambda_2 = 0.0695 \, \text{mm}^{-1}, \quad \lambda_3 = 0.0418 \, \text{mm}^{-1}, \\ \Lambda_1 &= 2.41 \, \text{mm}, \quad \Lambda_2 = 6.15 \, \text{mm}. \end{aligned} \qquad (15) $$

The corresponding granulometries are shown in Fig. (10). The values $\Lambda_1$ and $\Lambda_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.

<figure><img src="image_11.png" /><figcaption>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).</figcaption></figure>

### 4.2.3 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 $\eta > 0$ and $\vartheta > 0$:

$$ \gamma(x; \eta, \vartheta) = \frac{x^{\eta-1}}{\vartheta^{\eta} \Gamma(\eta)} e^{-x/\vartheta}, $$

where $\Gamma$ is the Gamma function, i.e. $\Gamma(t) = \int_{0}^{+\infty} dt \, t^{x-1}e^{-t}$. The expectation of the gamma law is $\eta\vartheta$. The distribution in number $f^X$ of the random variable $X = V/V_m$
was fitted by the gamma law (Marthinsen, 1996):

$$f_X(x) \approx \gamma \left( x; k, \frac{1}{k} \right) = \frac{x^{k-1} k^k}{\Gamma(k)} e^{-kx}, \qquad (16)$$

where $k = 5.56$. The granulometry in volume $g_X$ corresponding to $f_X$ is $g_X(x) = \varrho V f_X(x)$ where $\varrho > 0$ is chosen so that $g_X$ is of mean 1. Making use of $\int dx x f_X^X(x) = 1$, we find that $\varrho = 1/V_m$, and:

$$g_X(x) = x f_X(x) \approx \gamma \left( x; k + 1, \frac{1}{k} \right). \qquad (17)$$

Similarly, the distribution in volume of the random variable $V$ is:

$$g'_{X}(x) = \frac{1}{V_{m}} g_{X} \left( \frac{x}{V_{m}} \right) \approx \gamma \left( x; k + 1, \frac{V_{m}}{k} \right). \qquad (18)$$

For comparison with the concrete formulation, this distribution should be expressed in terms of radius $r$. We use:

$$G^V(r; V_m) \approx g'_X \left( \frac{4}{3} \pi r^3 \right), \qquad (19)$$

where we have assumed $V$ is a sphere. The expression in (13) is found to be minimal when:

$$\begin{aligned} V_m^1 &= 0.1157 \, \text{mm}^3, & V_m^2 &= 187.9 \, \text{mm}^3, & V_m^3 &= 2077 \, \text{mm}^3, \\ \Lambda_1 &= 2.01 \, \text{mm}, & \Lambda_2 &= 5.08 \, \text{mm}. \end{aligned} \qquad (20)$$

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.

## 5.1 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
<figure><img src="image_12.png" /><figcaption>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).</figcaption></figure>

its bounding box. This library, detailed in appendix (A), is generated with a Poisson intensity resulting in an average of $N = 200$ planes in the simulated volume. It is a full (non-truncated) Poisson granulometry. We wish to choose Poisson polyhedra following the truncated granulometry [$\Lambda_1, \Lambda_2$] with Poisson intensity $\lambda_I$ and volume fraction $f_I$. We first apply an homothety of factor $\alpha$ to all polyhedra to account for the change of scales between the library and that of the simulation:

$$R'_{i} = \alpha R_{i}, \qquad \alpha = \frac{N}{2\pi\sqrt{3}\lambda_{I}L}, \qquad (21)$$

where $L$ is the size of the simulated microstructure and $R_i$ is the polyhedron radius. The quantity $2\pi\sqrt{3}\lambda_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 $\Lambda_1 < R'_i < \Lambda_2$. To achieve the volume fraction $f_I$, we choose randomly in the retained Poisson polyhedra $n_p$ aggregates with:

$$n_p = \frac{f_I L^3}{K_0^T}, \qquad (22)$$

where $K_0^T$ is the average Poisson polyhedra volume in the truncated granulometry. In turn, we compute $K_0^T$ numerically using the formula of Miles-Lantuéjoul, to correct for border effects (see appendix A):

$$K_0^T = \alpha^3 \frac{1}{\sum_{i \in I} P[B^i \subset D]^{-1}} \sum_{i \in I} \frac{1}{P[B^i \subset D]} V^i, \qquad (23)$$

where $\alpha^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 \subset D]$ is given in appendix (A).

## 5.2 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 be-
long 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 cri-
terion 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** uni-
formly 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 $\varpi$, 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 trans-
lated in the direction **v** by steps $\varpi$ until there is a non empty intersection. When there
is intersection, the polyhedron is moved once in the opposite direction, with a transla-
tion length $\varpi$. 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 mech-
anism 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 vecto-
rially 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).

## 5.3 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 $\pi = 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.
<table><caption>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.</caption><tr><td>Class (mm)</td><td>6.15 – 12.5</td><td>2.4 – 6.15</td><td>0.05 – 2.4</td><td>0 – +∞</td></tr><tr><td>λ (mm<sup>-1</sup>)</td><td>0.0418</td><td>0.0695</td><td>0.518</td><td>0.518</td></tr><tr><td>Resolution (mm/voxel)</td><td>0.6</td><td>0.31</td><td>0.025</td><td>0.03</td></tr><tr><td>Reached filling ratio</td><td>0.38</td><td>0.40</td><td>0.42</td><td>0.43</td></tr><tr><td>Retained polyhedra</td><td>22.5%</td><td>41.9%</td><td>96.6%</td><td>100%</td></tr></table>
5.4 Microstructure generation

We generate a microstructure containing 1600³ voxels according to the parameters given in Tab. (2). The number of aggregates is about 10⁶. 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. The resulting overall aggregates volume fraction is f_t ≈ 65.2%, not far from the formulation (68.7%). 2D cross-sections of the microstructure are represented in Fig. (12).
<table><caption>Table 2: Parameters used for generating a three-scales concrete microstructure: Poisson intensity $\lambda_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 \approx 65.2\%$ and $f_t^{LPE} \approx 64.3\%$ after polyhedra disconnection.</caption><tr><td>Class (mm)</td><td>6.15 – 12.5</td><td>2.4 – 6.15</td><td>0.05 – 2.4</td></tr><tr><td>λ<sub>I</sub> (mm<sup>-1</sup>)</td><td>0.0418</td><td>0.0695</td><td>0.518</td></tr><tr><td>f<sub>I</sub></td><td>0.305</td><td>0.123</td><td>0.24</td></tr><tr><td>Resolution</td><td colspan="3">0.09 mm/voxel</td></tr><tr><td>Image size</td><td colspan="3">1600<sup>3</sup> voxels (144<sup>3</sup> mm<sup>3</sup>)</td></tr><tr><td>Step size ϖ</td><td colspan="3">2 voxels</td></tr></table>
<figure><img src="image_13.png" /><figcaption>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).</figcaption></figure>

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 \approx 65.2\%$ to $f_t^{LPE} \approx 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.

## 5.5 Model validation using covariance and granulometries measurements

We first check the influence of the step parameter $\bar{\omega}$ on the implantation algorithm with regards to covariances. We consider microstructures containing two classes of polyhedra with the largest sizes and varying step sizes $\bar{\omega} = 1, 2, 4$ voxels (Tab. 3). A 2D section of the microstructure with $\bar{\omega} = 2$ is represented in Fig. (13), and the covariances of each model are shown in Fig. (14). The three covariances are very
close to each other. Additionally, we verify that the aggregates percolate in the three microstructures, irrespective of the value of $\omega$. Accordingly, the step size parameter does not influence the covariance, and can not be used to control percolation.
<table><caption>Table 3: 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.).</caption><tr><td>Class (mm)</td><td>6.15 – 12.5</td><td>2.4 – 6.15</td></tr><tr><td>λ<sub>I</sub> (mm<sup>-1</sup>)</td><td>0.0418</td><td>0.0695</td></tr><tr><td>f<sub>I</sub></td><td>0.305</td><td>0.123</td></tr><tr><td>Resolution</td><td colspan="2">0.24 mm/voxel</td></tr><tr><td>Image size</td><td colspan="2">1024<sup>3</sup> voxels (246<sup>3</sup> mm<sup>3</sup>)</td></tr><tr><td>Step size ϖ</td><td colspan="2">1, 2, 4 voxels</td></tr></table>
<figure><img src="image_14.png" /><figcaption>Figure 13: 2D cross-section of the microstructure corresponding to the parameters given in Tab. (3) with step size $\omega = 2$.</figcaption></figure>

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.

# 6 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
<figure><img src="image_15.png" /><figcaption>Figure 14: Centered reduced covariance $C_{c,r}(h)$ of the three microstructures detailed in Tab. (3).</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>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).</figcaption></figure>

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 developed 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 optimizing 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).
<figure><img src="image_17.png" /><figcaption>Figure 16: Granulometry measurement on microstructures simulated according to the parameters given in Tab. (2) (blue) vs. theoretical granulometry (black).</figcaption></figure>

# 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 Algorithm implementation

### A.1.1 Space tessellation

A Poisson tessellation of a cubic subdomain $D$ of size $L$ included in the space $\mathbb{R}^3$ is first generated, following section (4.1). We limit ourselves to a Poisson points process in $S_3 \times [0; r_1]$, where $r_1 = \sqrt{3}L/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 $\lambda$ and to the measure of the set $S_3 \times [0; r_1]$:

$$N = 2\pi\sqrt{3}\lambda L. \qquad (24)$$

We implant $N_P$ planes, where $N_P$ follows a Poisson distribution (Eq. 10) with expectation $N$ above. The planes $(P_i)_{i \le 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, \sqrt{3}L/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 $\mathbf{x}/\|\mathbf{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. Additionally, only polyhedra that are not cut by the boundary of $D$ are labeled. The tessellation is defined as a set of planes $(\Pi_i)_{i \le N_p}$ 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 \le N_p}$ such that all points $(x, y, z)$ in the polyhedron verify:

$$s_i = \operatorname{sign}(a_i x + b_i y + c_i z + d_i), \qquad (25)$$

We now consider all straight lines $(\Delta_{i,j})_{i,j \le N_p}$ at the intersection of pairs of planes $\Pi_i$ and $\Pi_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 $\Delta_{i,j}$ among the remaining ones and its intersections with all planes in the tessellation $(\Pi_{k_c})_{k_c \le N_p}$. We select the $M-1$ intersecting points $(I_k)_{1 \le k \le 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 $\Delta_{i,j}$ with $D$, so that $\Delta_{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 \ne i, k \ne 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 $(\pm 1, \pm 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 whether 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 polyhedron 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} \leftarrow \min(x_{\min}, x_{I_k}, x_{I_{k+1}})$ and similarly for $x_{\max}$ and the other coordinates. Finally, the set of faces are updated by adding (if not already present) the planes $\Pi_i$ and $\Pi_j$ as well as the two planes that correspond to the intersection points $I_K$ and $I_{K+1}$.
<figure><img src="image_18.png" /><figcaption>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 $\Delta$, $\Delta$ 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.</figcaption></figure>

### A.1.3 Miles-Lantuéjoul correction

Discarding polyhedra that cross the boundary of the domain $D$ favor smaller polyhedra. The probability one polyhedron has to cross the domain boundary indeed decreases with its volume. To correct this bias, the Miles-Lantuéjoul correction consists in assigning to each polyhedron the weight $1/P[B \subset D]$, instead of equiprobability, where $P[B \in D]$ is the probability that the bounding box $B$ of an object is fully enclosed in $D$ (Serra, 1982). It is given by:

$$P[B \subset D] = \frac{(L - L_x)(L - L_y)(L - L_z)}{L^3}, \qquad (26)$$

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 shown in Fig. (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><img src="image_19.png" /><figcaption>Figure 18: Miles-Lantueéjoul correction: probability for an object to have its bounding box B entirely included in a cubic domain D.</figcaption></figure>

## A.2 Validation of the algorithm and polyhedra library

### A.2.1 Geometrical covariogram

We recall the expression for the geometrical covariogram $K(\mathbf{h})$ of a set $A$:

$$K(\mathbf{h}) = \mu_3(A \cap A_{-\mathbf{h}}), \qquad (27)$$

where $\mu_3$ is the Lebesgue measure on $\mathbb{R}^3$, and $A_{-\mathbf{h}}$ is the translation of $A$ by the vector $\mathbf{h}$, as shown in Fig. (21). For Poisson polyhedra of Poisson planes intensity $\lambda$:

$$K(h) = K(0)e^{-\pi\lambda\|\mathbf{h}\|}, \quad K(0) = \frac{6}{\pi^4\lambda^3}. \qquad (28)$$

where $K(0)$ is the volume expectation and $\|\mathbf{h}\|$ is the norm of $\mathbf{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 $\mathbf{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 \le N_t}$ is thus generated. The geometrical covariogram is estimated by averaging:

$$K(h) = \frac{1}{N_t} \sum_{i=1}^{N_t} \mu_3(A^i \cap A_{-h}^i). \qquad (29)$$

We also compute:

$$K^{\text{corr}}(h) = \frac{1}{\sum_{i} P[B^i \subset D]^{-1}} \sum_{i=1}^{N_t} \frac{1}{P[B^i \subset D]} \mu_3(A^i \cap A_{-h}^i). \qquad (30)$$
<figure><img src="image_20.png" /><figcaption>Figure 19: 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).</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>Figure 20: 3D visualizations of some Poisson polyhedra.</figcaption></figure>

where $B^i$ is the bounding box of the polyhedron $A^i$, as a corrected estimate of $K(h)$, following the Miles-Lantuéjoul correction (appendix A.1.3). Estimates of the geometrical covariograms $K(h)$ and $K^{corr}(h)$ are compared to the theoretical geometrical covariogram in Fig. (22) when $N = 200$. As expected, estimates are closer the theoretical values at smaller values of $h$. The corrected estimate also significantly improves the measure.

Relative errors using the corrected estimate for the average volume $K^{corr}(0)$ with respect to the theoretical value are given in Tab. (4) for increasing values of $N$. The average $\langle L/L_x \rangle$, which indicate the number of polyhedron along one dimension of $D$ is also given. At low values of $N$, the Poisson tessellation is not representative, resulting in large errors for the estimate of $K^{corr}(0)$. The value $N = 200$ offers a good compromise between computational costs and precision, with a relative error of 7%. For this value, $\langle L/L_x \rangle \approx 50$, which is consistent with measures in 2D with polygons, using the Miles-Lantuéjoul correction (C. Lantuéjoul, private communication, 2011).
<figure><img src="image_22.png" /><figcaption>Figure 21: Measure of the geometrical covariogram $K(\mathbf{h})$ as the volume of $A \cap A_{-\mathbf{h}}$, shown in gray for a given vector $\mathbf{h}$.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>Figure 22: Theoretical value and numerical estimates of the geometrical covariogram $K(h)$ and $K^{\text{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.</figcaption></figure>
<table><caption>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 $\langle L/L_x \rangle$ is given on row 3.</caption><tr><td>N</td><td>80</td><td>100</td><td>150</td><td>200</td><td>300</td></tr><tr><td>e<sub>rel</sub>(K<sup>corr</sup>(0))</td><td>73%</td><td>43%</td><td>14%</td><td>7%</td><td>6%</td></tr><tr><td>⟨L/L<sub>x</sub>⟩</td><td>19</td><td>23</td><td>35</td><td>47</td><td>70</td></tr></table>
### A.2.2 Covariance of a Poisson mosaic

We 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$ (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:

$$C(h) = fr(h) + f^2 [1 - r(h)], \qquad (31)$$

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:

$$C^*(h) = \frac{K_{A \cap D}(h)}{K_D(h)}, \qquad (32)$$

with $K_D(h) = L^2(L-h)$. As many as $N_A = 5 \times 10^4$ random points uniformly distributed in a domain $D$ are used to estimate $K_{A \cap D}(h)$ in vector images, without binarization. It is estimated by:

$$K_{A \cap D}(h) = \frac{1}{N_A} \sum_{i=1}^{N_A} g(\mathbf{x}_i)g(\mathbf{x}_i + \mathbf{h}), \qquad (33)$$

where $\mathbf{h}$ is a vector of norm $h$ oriented along the direction $x$, and $g(\mathbf{x}) = 1$ (resp. $g(\mathbf{x}) = 0$) if $\mathbf{x} \in A$ (resp. $\mathbf{x} \in A$). Additionally, 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><img src="image_24.png" /><figcaption>Figure 23: 2D cross-section of a 3D binary Poisson mosaic of size 256³ and parameters $N = 200$ and $f = 0.25$.</figcaption></figure>
<figure><img src="image_25.png" /><figcaption>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.</figcaption></figure>

### A.2.3 Poisson polyhedra library

We generate a library of $5 \times 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 $\lambda$ requires the rescaling given in section (5.1). This is computed analytically. For instance the plane $\Pi_i$ defined by $a_ix+b_iy+c_iz+d_i=0$ is implanted in point $(x_c, y_c, z_c)$ and rescaled by $\alpha$ by changing $d_i$ into:

$$d'_i = \alpha d_i - ax_c - by_c - cz_c + \frac{\alpha}{2} [a(x_{min} + x_{max}) + b(y_{min} + y_{max}) + c(z_{min} + z_{max})]. \quad (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 $\mathbf{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:

$$\max_{\mathbf{x}_c, R} \{R; \mathbf{V}_i^T \mathbf{x}_c + R \|V_i\| \le v_i, i = 1, \dots, m\}, \quad (35)$$

where $\mathbf{V}_i$ is a vector of dimension 3 such as $\mathbf{V}_i = (-s_ia_i, -s_ib_i, -s_ic_i)$ and with $v_i$ given by $v_i = s_id_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):

$$G(r) = 1 - \left( 1 + 4\pi\lambda r + \frac{\pi^4}{6}\lambda^2 r^2 + \frac{2}{9}\pi^5\lambda^3 r^3 \right) e^{-4\pi\lambda r} \quad (36)$$

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.

<figure><img src="image_26.png" /><figcaption>Figure 25: Cumulative granulometry of Poisson polyhedra: theoretical opening granulometry vs. numerical granulometry, computed by Chebyshev measures on the Poisson polyhedra library.</figcaption></figure>

### 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 $\theta$. 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):

$$f = 1 - e^{-\theta K_0}, \quad Q(h) = (1 - f)^{2-r(h)}, \quad (37)$$
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:

$$V^0 = \frac{1}{\sum_j P[B_j \subset D]^{-1}} \sum_j \frac{1}{P[B_j \subset D]} V_j, \quad (38)$$

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 = \alpha^3 V^0$, where $\alpha$ 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.

<figure><img src="image_27.png" /><figcaption>Figure 26: 2D cross-section of a Poisson polyhedra Boolean media with intensity $\lambda = 0.045$ and volume fraction $f = 0.2$. The image is discretized in a 500³ voxels grid.</figcaption></figure>

We measure the covariance of the complementary set of a Poisson polyhedra Boolean media with intensity $\lambda = 0.045$ and volume fraction $f = 0.2$. The 3D image is discretized on 500³ voxels. The measured covariance and theoretical expression are compared in Fig. (27). Excellent agreement is found between the two. For values of $\lambda$ smaller than $\lambda = 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.
<figure><img src="image_28.png" /><figcaption>Figure 27: Covariance $Q(h)$ of the complementary set of a Poisson polyhedra Boolean model with intensity $\lambda = 0.045$ and volume fraction $f = 0.2$: comparison between measurement and theoretical expression.</figcaption></figure>

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):

$$N_t = \frac{\pi N^3}{144\sqrt{3}}. \qquad (39)$$

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.

## References

* Bary, B., Ben Haha, M., Adam, E. & Montarnal, P. (2009) Numerical and Analytical Effective Elastic Properties of Degraded Cement Pastes. Cem. Concr. Res. 39, 902–912.

* Bazant, Z.B., Tabbara, M.R., Kazemi, M.T. & Pijaudier-Cabot, G. (1990) Random particle model for fracture of aggregate of fiber composites. J. Eng. Mech. 116, 1686–1705.

* Bentz, D.P. (1997) Three-dimensional computer simulation of Portland cement hydration and microstructure development. J. Am. Ceram. Soc. 80, 3–21.
Bentz, D.P., Garboczi, E.J., Haecker C.J. & Jensen, O.M. (1999) Effects of cement particle size distribution on performance properties of portland cement-based materials. *Cem. Concr. Res.* 29, 1663–1671.

Bernard, F., Kamali-Bernard S. & Prince, W. (2008) 3D multi-scale modelling of mechanical behaviour of sound and leached mortar. *Cem. Concr. Res.* 38, 449–458.

Beucher, S. & Meyer, F. (1992) *Mathematical Morphology in Image Processing*, Chap. The Morphological Approach to Segmentation: The Watershed Transformation, CRC Press (Marcel Dekker pub.), 433–481.

Bishnoi, S. & Scrivener, K.L. (2009) *µic*: A new platform for modelling the hydration of cements. *Cem. Concr. Res.* 39, 266–274.

Bullard, J.W. & Garboczi, E.J. (2006) A model investigation of the influence of particle shape on portland cement hydration. *Cem. Concr. Res.* 36, 1007–1015.

Comput. Struct. Caballero, A., López, C.M. & Caroll, I. (2006) 3D meso-structural analysis of concrete specimens under uniaxial tension. *Comput. Method Appl. M.* 195, 7182–7195.

Dequiedt, A.S., Coster, M., Chermant, J.L. & Jeulin, D. (2001) Towards a model of concrete mesostructure. *Cement Concrete Comp.* 23, 289–297.

Dunant, C.F., Bary, B., Giorla, A.B., Péniguel, C., Sanahuja, J., Toulemonde, C., Tran, A.B., Willot, F. & Yvonnet, J. (2013) A critical comparison of several numerical methods for computing effective properties of highly heterogeneous materials. *Adv. Eng. Softw.* 58, 1–12.

Escoda, J. (2012) *Modélisation morphologique et micromécanique 3D de matériaux cimentaires*. PhD Thesis, Mines ParisTech.

Escoda, J., Jeulin, D. & Willot, F. (2011) Simulation of 3D granular media by multiscale random polyhedra. In *Proceedings of the 13th International Congress for Stereology*, Guoquan Liu et al. (eds.), Beijing (electronic document, online at http://hal.archives-ouvertes.fr/hal-00879260).

Escoda, J., Willot, F., Jeulin, D., Sanahuja, J. & Toulemonde, C. (2011) Estimation of local stresses and elastic properties of a mortar sample by FFT computation of fields on a 3D image. *Cem. Concr. Res.* 41, 542–556.

Gal, E. & Kryvoruk, R. (2011) Meso-scale analysis of frc using a two-step homogenization approach. *Comput. Struct.* 89, 921–929.

Gallucci, E., Scrivener, K., Groso, A., Stampanoni, M. & Margaritondo, G. (2007) 3D experimental investigation of the microstructure of cement pastes using synchrotron X-ray microtomography (µCT). *Cem. Concr. Res.* 37, 360–368.
Garboczi, E.J. (2002) Three-dimensional mathematical analysis of particle shape using x-ray tomography and spherical harmonics: Application to aggregates used in concrete. *Cem. Concr. Res.* 32, 1621–1638.

Garboczi, E.J. & Bentz, D.P. (1992) Computer simulation of the diffusivity of cement-based materials. *J. Mater. Sci.* 27, 2083–2092.

Gilbert, E.N. (1962) Random subdivisions of space into crystals. *Ann. Math. Stat.* 33, 958–972.

Grant, M. & Boyd, S. (2008) Graph implementations for nonsmooth convex programs. In V. Blondel, S. Boyd and H. Kimura, editors, *Recent advances in Learning and Control*, Lecture Notes in Control and Information Sciences 371, 95–110. Springer-Verlag Limited.

Grant, M. & Boyd, S. (2011) CVX: Matlab software for disciplined convex programming, version 1.21. http://cvxr.com.

Haecker, C.J., Garboczi, E.J., Bullard, J.W., Bohn, R.B., Sun, Z., Shah, S.P. & Voigt, T. (2005) Modeling the linear elastic properties of Portland cement paste. *Cem. Concr. Res.* 35, 1948–1960.

Häfner, S., Eckardt, S., Luther, T. & Könke, C. (2006) Mesoscale modeling of concrete: Geometry and numerics. *Comput. Struct.* 84, 450–461.

Hain, M. & Wriggers, P. (2007) Numerical homogenization of hardened cement paste. *Comput. Mech.* 42, 197–212.

He, H., Guo, Z., Stroeven, P., Stroeven M. & Sluys, L.J. (2010) Strategy on simulation of arbitrary-shaped cement grains in concrete. *IAS* 29, 79–84.

Jeulin, D. (1991) Modèles morphologiques de structures aléatoires et de changement d'échelle. Thèse de doctorat d'état, Université de Caen.

Kapur, J.N., Sahoo, P.K. & Wong, A.C.K. (1985) A new method for gray-level picture thresholding using the entropy of the histogram. *Graph. Model. Im. Proc.* 29, 273–285.

Kim, S.M. & Abu Al-Rub, R.K. (2011) Meso-scale computational modeling of the plastic-damage response of cementitious composites. *Cem. Concr. Res.* 41, 339–358.

Knuth, D.E. (1997) *The art of computer programming, vol. 2: Seminumerical algorithms*, Addison-Wesley (3rd edition), Boston.

Kwan, A.K.H., Wang, Z.M. & Chan, H.C. (1999) Mesoscopic study of concrete ii: nonlinear finite element analysis. *Comput. Struct.* 70, 545–556.

Marthinsen, K. (1996) Comparative analysis of the size distributions of linear, planar, and spatial Voronoi cells. *Mater. Charact.* 36, 53–63.
Matheron, G. (1971) Les polyèdres Poissoniens isotropes. In *Actes de conférence du troisième colloque Européen sur la fragmentation* 2 (Dechema ed.), 509–534, Cannes.

Matheron, G. (1972) Ensembles fermés aléatoires, ensembles semi-markoviens, et polyèdres poissoniens. *Adv. Appl. Probab.* 4, 508–541.

Matheron, G. (1975) *Random sets and integral geometry*, Wiley, New-York.

Nagai, G., Yamada, T. & Wada, A. (1998) Stress analysis of concrete material based on geometrically accurate finite element modeling. *Third International Conference on Fracture Mechanics of Concrete and Concrete Structures (FRAMCOS-3)*, AEDI-FICATIO publishers, P 42, 1077–1086.

Nelder, J.A & Mead, R. (1965) A simplex method for function minimization. *Comput. J.* 7, 308–313.

Otsu, N. (1979) A threshold selection method from gray-level histograms. *IEEE Trans. Syst. Man. Cybern.* 9, 62–66.

Pedersen, R.R., Simone, A., Stroeven, M. & Sluys, L.J. (2007) Mesoscopic modelling of concrete under impact. In *Proceedings of the 6th International Conference on Fracture Mechanics of Concrete and Concrete Structures*, A. Carpinteri, P. G. Gambarova, G. Ferro and G. Plizzari (eds.) Catania, Italy.

Quenec'H, J.L., Coster, M., Chermant, J.L. & Jeulin, D. (1992) Study of the liquid-phase sintering process by probabilistic models: Application to the coarsening of WC-Co cermets. *J. Microsc.* 168, 3–14.

Schlangen, E. & van Mier, J. (1992) Simple lattice model for numerical simulation of fracture of concrete materials and structures. *Mater. Struct.* 25, 534–542.

De Schutter, G. & Taerwe, L. (1993) Random particle model for concrete based on Delaunay triangulation. *Mater. Struct.* 26, 67–73.

Serra, J. (1982) *Image Analysis and Mathematical Morphology*, vol. 1, Academic Press (publ.), New-York.

Sezgin, M. & Sankur, B. (2004) Survey over image thresholding techniques and quantitative performance evaluation. *J. Electron. Imaging* 13, 146–165.

Stark, J. (2011) Recent advances in the field of cement hydration and microstructure analysis. *Cem. Concr. Res.* 41, 666–678.

Stroeven, P. (2000) A stereological approach to roughness of fracture surfaces and tortuosity of transport paths in concrete. *Cement Concrete Comp.* 22, 331–341.

Stroeven, P., Le, L.B.N., Sluys, L.J. & He, H. (2012a) Porosimetry by double-random multiple tree structuring in virtual concrete. *IAS* 31, 55–63.
Stroeven, P., Le, L.B.N., Sluys, L.J., & He, H. (2012b) Porosimetry by random node structuring in virtual concrete. *IAS* 31, 79–87.

Thomas, J.J., Biernacki, J.J., Bullard, J.W., Bishnoi, S., Dolado, J.S., Scherer, G.W. & Luttge, A. (2011) Modeling and simulation of cement hydration kinetics and microstructure development. *Cem. Concr. Res.* 41, 1257–1278.

Šmilauer, V. & Bažant, Z. (2010) Identification of viscoelastic C-S-H behavior in mature cement paste by FFT-based homogenization method. *Cem. Concr. Res.* 40, 197–207.

Wang, Z.M., Kwan, A.K.H. & Chan, H.C. (1999) Mesoscopic study of concrete I: generation of random aggregate structure and finite element mesh. *Comput. Struct.* 70, 533–544.

Willot, F., Abdallah, B. & Pellegrini, Y.P. (2014) Fourier-based schemes with modified Green operator for computing the electrical response of heterogeneous media with accurate local fields. *Int. J. Numer. Meth. Eng.* 98, 518–533.

Willot, F. & Jeulin, D. (2009) Elastic behavior of composites containing boolean random sets of inhomogeneities. *Int. J. Eng. Sci.* 47, 313–324.

Wriggers, P. & Moftah, S.O. (2006) Mesoscale models for concrete: Homogenisation and damage behaviour. *Finite Elem. Anal. Des.* 42, 623–636.

Ye, G., van Breugel, K. & Fraaij, A.L.A. (2003) Three-dimensional microstructure analysis of numerically simulated cementitious materials. *Cem. Concr. Res.* 33, 215–222.

Zaitsev, Y. & Wittmann, F. (1981) Simulation of crack propagation and failure of concrete. *Mater. Struct.* 14, 357–365.