# Representative volume elements for matrix-inclusion composites - a computational study on periodizing the ensemble

Matti Schneider, Marc Josien, Felix Otto

► To cite this version:

Matti Schneider, Marc Josien, Felix Otto. Representative volume elements for matrix-inclusion composites - a computational study on periodizing the ensemble. Journal of the Mechanics and Physics of Solids, 2022, 158, pp.104652. 10.1016/j.jmps.2021.104652. cea-03713963

HAL Id: cea-03713963

https://cea.hal.science/cea-03713963

Submitted on 5 Jul 2022

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

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

# Representative volume elements for matrix-inclusion composites - a computational study on the effects of an improper treatment of particles intersecting the boundary and the benefits of periodizing the ensemble

Matti Schneider<sup>a, *</sup>, Marc Josien<sup>b</sup>, Felix Otto<sup>c</sup>

## ARTICLE INFO

**Keywords:**

- Quantitative stochastic homogenization
- FFT-based computational homogenization
- Thermal conductivity
- Screening the boundary-layer error
- Microstructure modeling

## ABSTRACT

We investigate volume-element sampling strategies for the stochastic homogenization of particle-reinforced composites and show, via computational experiments, that an improper treatment of particles intersecting the boundary of the computational cell may affect the accuracy of the computed effective properties. Motivated by recent results on a superior convergence rate of the systematic error for periodized ensembles compared to taking snapshots of ensembles, we conduct computational experiments for microstructures with circular, spherical and cylindrical inclusions and monitor the systematic errors in the effective thermal conductivity for snapshots of ensembles compared to working with microstructures sampled from periodized ensembles.

We observe that the standard deviation of the apparent properties computed on microstructures sampled from the periodized ensembles decays at the scaling expected from the central limit theorem. In contrast, the standard deviation for the snapshot ensembles shows an inferior decay rate at high filler content. The latter effect is caused by additional long-range correlations that necessarily appear in particle-reinforced composites at high, industrially relevant, volume fractions. Periodized ensembles, however, appear to be less affected by these correlations.

Our findings provide guidelines for working with digital volume images of material microstructures and the design of representative volume elements for computational homogenization.

## 1. Introduction

### 1.1. State of the art

The concept of representative volume element (RVE) plays a central role for predicting the effective properties of random heterogeneous materials. Originally, Hill (1963) defined an RVE as a bounded domain that is statistically typical of the mixture and sufficiently large to render the effects of artificially imposed boundary conditions negligible. Drugan and Willis (1996) improved upon the original concept by relaxing the previous definition, that requires the RVE to be typical for all statistics of interest. Rather,

<footnote><sup>a</sup> Karlsruhe Institute of Technology (KIT), Institute of Engineering Mechanics, Germany</footnote>

<footnote><sup>b</sup> CEA, DES, IRESNE, DEC, Cadarache, F-13108, Saint-Paul-Lez-Durance, France</footnote>

<footnote><sup>c</sup> Max-Planck Institute for Mathematics in the Sciences (MPI MiS), Germany</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: matti.schneider@kit.edu (M. Schneider).</footnote>
they focused on the material properties and defined an RVE as a bounded domain whose apparent properties are sufficiently close to the effective properties of the infinitely large medium. This shift in paradigm, which introduced the infinite-volume limit as a suitable reference quantity, enabled much smaller RVEs to be used, with obvious advantages in terms of computational complexity. With this practical definition at hand, representativity of a volume element could be determined by investigating the effective properties computed on a sequence of volume elements with increasing size, see Gusev (1997) or Segurado and Llorca (2002). Kanit et al. (2003) enriched this approach by statistical ideas, proposing to monitor the apparent properties of several realizations of the random material on cells of fixed size. They provided theoretical arguments as well as empirical data in support of using several realizations of the ensemble on smaller volumes to arrive at the same accuracy as for a single “large” RVE. However, they also showed that care has to be taken for “small” cells, as there may be a bias caused by deterministic size effects which cannot be eliminated through ensemble averaging. Still, Kanit et al. assume the introduced bias to be small compared to the dispersion, i.e., the variance, of the effective properties on cells of fixed size. They propose to numerically identify the decay behavior of the dispersion, and to define an RVE as a (finite collection of) volume element(s) with suitably small dispersion.

From a practical point of view, it is imperative to minimize the computational effort associated to numerically determining the apparent properties of volume elements. Finding an RVE may be regarded as selecting a suitable termination criterion for an iterative procedure which either processes more realizations with the same cell size or increases the considered volume element size. A number of different approaches was reported in the literature. In the spirit of Kanit et al. (2003), one may aim at a satisfactory confidence in average quantities (Gitman et al., 2007; Pelissou et al., 2009; Ghossein and Lévesque, 2012; Harper et al., 2012), for instance by employing Huet’s partition theorem (Huet, 1990) to deduce confidence intervals of the apparent properties by statistical sampling. In addition to the effective coefficients (Salmi et al., 2012; Trias et al., 2006; Terada et al., 2000; Shan and Gokhale, 2002), effective energies (Saroukhani et al., 2015) or mean values and concentrations of the microstructural fields (Trias et al., 2006; Shan and Gokhale, 2002; Stroeven et al., 2004; Gitman et al., 2007, 2006) may serve as suitable average quantities. Alternatively, bounds of the effective properties, inferred from Dirichlet and Neumann boundary conditions or analytical estimates, may be monitored, and RVEs arise provided these bounds are sufficiently close (Ostoja-Starzewski, 2006; Salmi et al., 2012; Hoang et al., 2016). This strategy was also applied for inelastic (Jiang et al., 2001) materials and at finite strains (Khisaeva and Ostoja-Starzewski, 2006; Ma et al., 2015). Some authors include further statistical measures, like correlation functions (Trias et al., 2006; Shan and Gokhale, 2002; Niezgoda et al., 2010; Teferra and Graham-Brady, 2018) or deviations from isotropy (for isotropic media) (Salmi et al., 2012a; Moussaddy et al., 2013), into their RVE definition. Relying on a single criterion may lead to a premature termination (Moussaddy et al., 2013), and more recent works combine different termination criteria (Salmi et al., 2012; Moussaddy et al., 2013; Trias et al., 2006; Saroukhani et al., 2015) or rely upon specific adaptive strategies (Doškář et al., 2018) based on microstructural building blocks, the Wang tiles (Šedlbauer and Lepš, 2019).

The inferred RVE size depends on a number of factors, including the underlying physics, the microstructure geometry and the material contrast (Stroeven et al., 2004). In the case of high or even infinite material contrast (Dirrenberger et al., 2014) or non-linear behavior (Stroeven et al., 2004; Gitman et al., 2007; Pelissou et al., 2009), the influence of the microstructure geometry gets more pronounced, leading to larger RVE sizes or even precluding the existence of an RVE (Gitman et al., 2007). Consequently, any recommendation for the RVE size depends strongly on the considered material, and cannot be directly transferred to other materials (Matouš et al., 2017). In particular, studies on the RVE size need to be carried out for each material under consideration, adding importance to computationally efficient strategies for identifying RVEs. Numerous RVE studies were conducted for various material classes, including granular materials (Evesque, 2005), heterogeneous polymers (Jean et al., 2011; Figliuzzi et al., 2016; Mirkhalaf et al., 2016), concrete (Sebsadji and Chouicha, 2012; Escoda et al., 2016), particle-reinforced composites (El Moumen et al., 2015; Ghossein and Lévesque, 2015; Savvas et al., 2016), polycrystalline materials (Ranganathan and Ostoja-Starzewski, 2008; Qidwai et al., 2012; Bouchedjra et al., 2018), porous metals (Fritzen et al., 2012; Masson et al., 2015; Wang et al., 2015), random composites (Chen et al., 2019; El Moumen et al., 2021), electrode and anode materials (Prill et al., 2017; Abdallah et al., 2016), explosives (Gasnier et al., 2015) and materials used in the food industry (Kanit et al., 2006, 2011). Improvements in computational power and algorithmic advances gave significant impetus to modern multiscale methods in engineering, also dealing with strain-gradient materials (Kouznetsova et al., 2002, 2004) and materials with softening effects (Nguyen et al., 2010) or debonding (Swaminathan and Ghosh, 2006), see Matouš et al. (2017) for a recent overview. In general, one may distinguish two succinctly different strategies for acquiring the microstructures necessary for multiscale methods in engineering. On the one hand, synthetic microstructures may be generated based on statistical data available on the material under consideration, see Bargmann et al. (2018) for a recent overview. On the other hand, digital volume images of real microstructures (Landis and Keane, 2010; Guven and Cinar, 2019) may be obtained directly, serving as what we call snapshots of the random material. Both approaches complement each other, and may be used cooperatively. Indeed, synthetic microstructures offer full control over the ensemble, and thorough RVE studies can be conducted. However, the relation of these synthetic structures to their real counterparts is not always apparent. In contrast, snapshots are real, but their acquisition is typically costly, which makes it hard to study their representativity and to quantify the underlying uncertainty.

Parallel to developments in engineering, homogenization theory was established as a rigorous foundation for the up-scaling of random heterogeneous media. Building upon previous work in the periodic setting (De Giorgi and Spagnolo, 1973; Babuska, 1973; Larsen, 1975), Papanicolaou and Varadhan (1981) and Kozlov (1978) established qualitative results for stochastic homogenization for stationary and ergodic ensembles of coefficient fields. In this context and informally speaking, stationarity means that the statistical properties of the ensembles are translation-invariant, whereas ergodicity of a random function means that the variance of its average over domains with increasing volume vanishes as this volume becomes infinite (Lantuéjoul, 1991), ensuring the convergence of space averages towards expectations for samples with increasing sizes. In particular, as a consequence of Birkhoff's
theorem (Zhikov et al., 1994, Thm. 7.2), ergodicity implies that the effective properties may be determined from a single, infinitely large realization of the microstructure (almost surely).

Classically, the effective properties in stochastic homogenization are defined by spatially averaging so-called corrector fields which are defined on the whole space. The RVE methodology, i.e., working with cells of large but finite size, is imperative from a computational point of view. The error of the apparent property associated to a single finite-sized cell compared to the true effective properties of stochastic homogenization naturally decomposes into two parts (Gloria et al., 2015, Eq. (13)). The random error quantifies the fluctuation of the apparent properties on cells of fixed size, whereas the systematic error measures the deviation of the expectation of the apparent property on cells of fixed size from the true effective property. These two error contributions correspond to the bias and dispersion considered by Kanit et al. (2003). Theoretical guarantees that the systematic as well as the random error converge to zero as the cell size goes to infinity were provided by Sab (1992), Bourgeat and Piatnitski (2004) and Owhadi (2003) for Dirichlet, Neumann and periodic boundary conditions imposed on the boundary of the considered cells.

These results hold in the general setting of stationary and ergodic random media. Thus, they may only provide qualitative convergence results, i.e., convergence of the random and the systematic error to zero, without rates. To obtain specific convergence rates, the correlation properties of the coefficient field need to be quantified, see Section 2 below for a discussion. Theoretical works also shed some light on the two different approaches for obtaining microstructures discussed earlier. For the mathematical analysis, the full ensemble of microstructures needs to be given, and its stochastic properties are encoded by the underlying correlation structure. Acquiring snapshots corresponds to drawing a specific sample of the ensemble, and restricting the microstructure to a cell of finite size. This approach is straightforward, and is tacitly assumed in the qualitative convergence results (Sab, 1992; Bourgeat and Piatnitski, 2004; Owhadi, 2003). Generating synthetic microstructures is more subtle, as those are typically generated on a rectangular cell with periodic boundary conditions, i.e., on a torus, mathematically speaking. Unlike the ensembles used in the snapshot strategy, these periodized ensembles retain the stationarity property of the whole-space ensemble. However, ergodicity is lost. The RVE methodology builds upon the fact that the loss of ergodicity fades away as the size of the volume element goes to infinity. Synthetically generated elements are more than simple restrictions of the infinite-space ensemble, but require the underlying “construction plan” to be known (Egloffe et al., 2015, Sec. 3.2.3), for instance via pushing forward a specific probability distribution by a (deterministic) mapping. With this more fundamental knowledge at hand, it is natural to draw samples from suitably periodized ensembles for the finitely sized cells, first formalized by (Sab and Nedjar, 2005). In practical engineering terms, microstructure-generation algorithms are - in principle - deterministic, but some of their input parameters are drawn from probability distributions. Combining the random input and the deterministic algorithm gives rise to the ensemble, and drawing samples from the periodized ensemble corresponds to working on the torus instead of a box in Euclidean space, for instance when computing distances or applying filters. From a theoretical point of view, the periodized ensemble arises by drawing samples from the whole-space ensemble conditioned on periodicity of the realizations (with previously fixed period). This interpretation, however, is of little use for practical algorithms.

Over the years, in terms of theoretical reasoning as well as numerical evidence, it has been realized that the snapshot strategy may introduce non-negligible boundary effects compared to the periodization approach. This observation surfaced in numerical multiscale methods (Hou and Wu, 1997; Ming and Zhang, 2005), where the finite-element ansatz functions were modified in order to account for heterogeneity on a sub-element scale by solving a corrector-type equation with suitable boundary conditions. In the setting of periodic homogenization it was realized that working on cubic volume elements with an edge length $L$ that is not an integer multiple of the period leads to an error in the effective properties that is of order $1/L$. To reduce this so-called resonance error, oversampling/filtering techniques (Hou and Wu, 1997; Blanc and Le Bris, 2010) were introduced. Alternative approaches screen the boundary effect by modifying the cell problem, for instance by adding a “massive”, zero-th order term to the elliptic operator defining the corrector problem (Gloria, 2011, Thm. 1). In a similar spirit, screening strategies based on parabolic (Abdulle et al., 2019; Mourrat, 2019) or hyperbolic (Arjmand and Runborg, 2016) versions of the corrector problem have been considered.

In the context of stochastic homogenization, the snapshot approach also gives rise to a boundary-layer error, just as for periodic homogenization with non-matching period. In the stochastic setting, Egloffe et al. (2015, Eq. (3.4)) expected the systematic error to decay with the rate $1/L$, as for the periodic case. Recently, Clozeau et al. (2021) provided a theoretical argument for the $1/L$-scaling, at least in the asymptotic setting of vanishing material contrast. Because of the slow decay of the random boundary layer, filtering techniques are expected to have no effect on the scaling $1/L$ (Gloria, 2008, Sec. 4). However, screening by a massive term, in conjunction with extrapolation in the massive parameter, has been shown to reduce the systematic error to $1/L^d$ (Gloria et al., 2015, Thm. 2). Based on screening and extrapolation, Mourrat (Mourrat, 2019, Prop. 1.1 & Th. 1.2) devised a numerical algorithm that, on the basis of a snapshot, extracts the effective behavior up to the optimal total-error rate $1/L^{d/2}$ at the cost of on the order of $L^d$ operations. Hence, the most effective improvements in the decay of the boundary-layer error rest upon a modification of the corrector problem. However, from an engineering point of view, modifying the underlying partial differential equation for more general problems to be homogenized, e.g. elasticity or inelasticity, can only be a last resort, even if it improves the convergence rate. Apart from physical reasons, such an approach requires modifying existing computational homogenization codes at their core.

Recently, Clozeau et al. (2021) worked out the details of an attractive alternative to modifying the corrector problem. Instead, they modify the ensemble by periodizing it on cells of finite size, thus recovering the more favorable $1/L^d$ convergence rate of the systematic error (without modifying the corrector equation). This is in line with the existing bounds of Gloria et al. (2015) and Egloffe et al. (2015, Sec. 3.2.3), which concern i.i.d. random conductivities on a discrete lattice, and with Khoromskia et al. (2020), who investigated an ensemble based on a discrete Poisson point process numerically. However, periodizing the aforementioned ensembles is trivial; on the contrary, Clozeau et al. (2021) deal with a less academic and more realistic class of ensembles, namely
<figure><img src="image_4.png" /><figcaption>**Fig. 1.** Illustrating the disturbance in statistical homogeneity created by using periodic boundary conditions for a non-periodic snapshot (left) and by appropriately reflecting a periodic microstructure (right) for the geometries shown in Fig. 2.</figcaption></figure>

coefficient fields generated by applying a nonlinear transformation to a Gaussian field with (essentially) integrable correlation function. Moreover, Clozeau et al. (2021) do not only provide an estimate of the systematic error but actually derive the leading-order term in $L^{-d}$, which is generically non-vanishing.

For applications in engineering, the result (Clozeau et al., 2021) suggests that working with periodized ensembles leads to a favorable systematic error compared to the snapshot approach. This is very much in line with observations in engineering. For instance, considering Poisson–Voronoi microstructures, Kanit et al. (Kanit et al., 2003, Sec. 7) observed that "the bias [systematic error] introduced by the periodic boundary conditions is found to be much smaller than for the other boundary conditions". This observation was confirmed by other authors (Sab and Nedjar, 2005; Yue and W.E., 2007). For Dirichlet and Neumann boundary conditions, Salmi et al. (2012) reported faster convergence when using volume elements whose faces do not intersect with heterogeneities than for a volume elements which intersect the heterogeneities. In their paper on periodization, Sab and Nedjar (2005) noted that "the use of homogeneous strain or stress boundary conditions introduces a systematic bias in the estimation of the overall elasticity tensor because the corresponding $\Gamma$ operator does not preserve the statistical homogeneity of the medium". Notice the careful wording used by Sab and Nedjar (2005) - they deliberately distinguish periodic boundary conditions and the periodized ensemble. Indeed, it is possible to use Dirichlet or Neumann boundary conditions for a periodized ensemble (Ostoja-Starzewski and Schulte, 1996, 1998) and to use periodic boundary conditions for a non-periodic medium (Michel et al., 1999). Both approaches interfere with the statistical homogeneity of the medium. This is illustrated in Fig. 1(a), which shows the periodic extension of a two-dimensional non-periodic snapshot of a material with circular inclusions. Due to a mismatch at the boundary, additional non-circular particle shapes are introduced. Something similar happens when working with Dirichlet and Neumann boundary conditions for a periodic medium. To see this, let us consider the case of Neumann boundary conditions in isotropic conductivity. To compute the effective flux corresponding to a prescribed temperature gradient $\tilde{\xi}$, it is possible to use periodic boundary conditions on a doubled cell (with $\tilde{\xi}$ serving as one of the edges), see Fig. 1 for $\tilde{\xi} = e_x$, where the coefficients are reflected in the directions orthogonal to $\tilde{\xi}$ and periodically extended in $\tilde{\xi}$-direction. Similar to the periodization of a non-periodic microstructure, this reflection procedure introduces non-circular particles which disturb statistical homogeneity.

## 1.2. Contributions

This article presents quantitative results on the consequences of improper treatment of particles intersecting the boundary of computational volume elements when homogenizing particle-reinforced composites. From an interdisciplinary perspective, the study provides a quantification for observations made in the engineering community (Kanit et al., 2003; Sab and Nedjar, 2005; Salmi et al., 2012) and extends related studies in the mathematics community (Khoromskia et al., 2020; Khoromskia and Khoromskij, 2020) to microstructures of industrial relevance and complexity, utilizing real physical parameters, investigating spherical and cylindrical fillers. In particular, we do not claim to introduce a new method, but inspect the role of the microstructure-acquisition process within RVE-size studies closely.
<figure><img src="image_5.png" /><figcaption>(a) Taking a snapshot introduces cut particles (green) (b) For a properly periodized ensemble, all inclusions in addition to original particles (red) belong to the same species</figcaption></figure>

<caption>Fig. 2. Illustrating snapshot coefficients $A_L^{sn}$ vs. periodized coefficients $A_L^{per}$.</caption>

More precisely, to assess the effects of taking snapshots, we rely upon well-established microstructure-generation methods (Williams and Philipse, 2003; Schneider, 2017) for composites with spherical and cylindrical fillers (see Section 3 for details), which automatically generate periodic ensembles and take snapshots from sufficiently large periodic structures to emulate snapshots of real composites, Fig. 2. As discussed in Section 1.1, different criteria may be used for assessing the representativity of a volume element (or a collection thereof). To minimize the ambiguities and to ensure comparability to the estimates established in the mathematical community, we monitor the systematic and the random error for the two considered microstructure-synthesis protocols, see Section 2 for precise definitions and a streamlined presentation of quantitative stochastic homogenization. Note that we are primarily interested in computing the effective properties of the random material in an efficient way, and do not focus on additional statistics.

Both the systematic and the random error of the apparent properties require the expectation of the apparent properties on cells of fixed, but finite, size to be computed to high accuracy. Previous studies (Khoromskia et al., 2020; Khoromskia and Khoromskij, 2020) revealed that 10000 realizations were sufficient to extract the desired convergence rates. At this point, we would like to reiterate our word of warning that the study at hand should not be mistaken as a proposal for a practical method (in view of the 10000 samples). Rather, we would like to present an accurate study under “laboratory conditions”, quantifying the effects of boundary-particle treatment in particle-reinforced composites. For this purpose, we selected thermal conductivity as the physical model under consideration. The (realistic) material properties used have a low material contrast of six, and we utilize modern computational homogenization techniques based on the fast Fourier transform (FFT), pioneered by Moulinec and Suquet (1994, 1998). Some of our investigations concern a similar setup as Altendorf et al. (2014), but have a different focus. Whereas Altendorf et al. (2014) identify RVEs for short-fiber reinforced composites to engineering accuracy, balancing accuracy and efficiency, our investigation is focused on a very precise quantification of the uncertainty inherent to the periodized and the snapshot protocols. Note that we impose periodic boundary conditions for the snapshots as well, see Fig. 1 and the associated discussion.

After exposing the microstructure-generation tools in Section 3, the computational results are presented in Section 4, comprising over 250000 simulation runs, some with up to more than one billion degrees of freedom. We find that the systematic error for the periodized protocol is typically lower than for the snapshot approach. The random error shows a similar trend. Moreover, the decay behavior of the random error in the snapshot protocol differs from the periodized setting. We discuss the latter finding in Section 5 more thoroughly, identifying this behavior as a rather natural consequence of the considered setting.

While the ensembles considered here do not fall into the class covered by the theory – they neither have finite range<fnref n="1" /> of dependence nor a suitable functional inequality is known, see Duerinckx and Gloria (2020b,a) – we expect that they belong to the same universality class. Thus, the composites under consideration serve as benchmarks for the established quantitative homogenization theory, and are expected to stimulate further theoretical research.

<footnote n="1">¹ Finite range of dependence should not be confused with the notion of integral range from geostatistics (Matheron, 1989; Lantuéjoul, 2002), as the latter is weaker.</footnote>
From a practical engineering perspective, these results shed light on modern up-scaling techniques based on digital volume images, for instance obtained by micro-computed tomography. Indeed, the latter images naturally correspond to snapshots of ensembles with “improperly” treated particles at the boundaries. For the specified scenario (thermal conductivity with low contrast) and the large database at hand, the empirical probability of a computation on a single cell (of predefined size) to provide results sufficiently close to the effective properties may be deduced. We observe that, for this seemingly innocent setup, two accurate digits is perfectly reasonable for the snapshot protocol, whereas obtaining three significant digits requires rather large cells.

Last but not least, let us emphasize that periodization of ensembles does not represent a panacea curing the problems of computational homogenization. Indeed, there are ensembles which do not admit a periodization (Dirrenberger et al., 2014; Jeulin, 2016; Sukiman et al., 2017).

## 2. Theoretical background

We are concerned with a stationary and ergodic random field of thermal conductivity tensors A in d spatial dimensions, i.e., every realization is a field on $R^d$ with values in symmetric positive definite $d \times d$-matrices. In stochastic homogenization (Papanicolaou and Varadhan, 1981; Kozlov, 1978), the effective conductivity tensor $\bar{A}$ is sought, and may be determined for any prescribed (negative) temperature gradient $\bar{\xi} \in R^d$ via the expectation (or ensemble average) $\langle \cdot \rangle$ of the heat flux

$$ \bar{A}_{\bar{\xi}} = \langle A_{\bar{\xi}} \rangle, \qquad (2.1) $$

where $\bar{\xi}$ denotes the corresponding local (negative) temperature gradient, i.e.,

$$ \bar{\xi} = \bar{\bar{\xi}} + \Gamma \bar{\xi} \qquad (2.2) $$

holds in terms of the Helmholtz projector $\Gamma$ onto stationary gradient fields, and the associated heat flux $q = A\bar{\xi}$ is divergence-free, i.e., the equation

$$ \Gamma q = 0 \qquad (2.3) $$

is satisfied. Due to the ergodicity assumption, the effective conductivity tensor $\bar{A}$ may also be computed in terms of a single realization involving an infinite-volume limit, see Chapter 7 in Zhikov et al. (1994) for details.

For computational purposes it is necessary to work on finite domains. For a fixed cube $Q_L = [-\frac{L}{2}, \frac{L}{2}]^d$, assume that we are given an ensemble $A_L$ that generates coefficient fields on the cell $Q_L$. The latter can be obtained as snapshots of the random conductivity tensor field A, which is given on the whole space, on the cube $Q_L$. In this case, we will write $A_L^{sn}$ for the corresponding ensemble. An alternative approach proceeds by periodizing the ensemble (Sab and Nedjar, 2005). In this case, we will denote the ensemble by $A_L^{per}$. Taking snapshots of an ensemble is straightforward because it can be performed on given realizations, whereas periodization is a more subtle process (to be discussed below) and requires working on the level of the ensemble itself. However, periodization is quite natural when generating synthetic microstructures. In any case, we assume the periodized ensemble to be stationary (on the torus).

Let us illustrate what we mean by periodizing the ensemble. For i.i.d. discrete coefficient fields (Gloria et al., 2015) and the Poissonian ensemble (Gloria and Otto, 2016), this procedure is trivial. For Gaussian fields, periodization amounts to periodizing the correlation function, see Clozeau et al. (2021). For the article at hand, we consider a particle-reinforced composite with a single species of particles, for instance spherical particles, and specific macroscopic statistics. When taking a snapshot of the coefficient field on a box $Q_L$ (thus building $A_L^{sn}$), particles intersecting the boundary of $Q_L$ are cut, giving rise to another species of particles with a different shape, see Fig. 2(a). In particular, taking snapshots of the particles does not preserve the morphological statistics of the whole-space ensemble, in general.

To understand how periodization (Sab and Nedjar, 2005) proceeds let us look at the way an ensemble of particle-filled composites is built. In real life, the fillers do not overlap (unless this is part of the production process, for instance in sintering). Thus, when generating an ensemble, an overlap check needs to be made. For the periodized ensemble, the particles are placed in a fixed cell $Q_L$, and the overlap checks are made not only between the original particles, but also with their periodic translations, see Fig. 3 for an illustration. Following such a protocol with statistical properties properly matching the whole-space counterpart, see Fig. 2(b), gives rise to the periodized ensemble and the associated coefficient field $A_L^{per}$. In particular, due to the absence of the cutting process, no alteration of the particle shape takes place. As will be shown below, the cut particles of $A_L^{sn}$ are considered responsible for deteriorating the quality of the computed homogenized coefficient with respect to the periodized coefficient field $A_L^{per}$ (see Eq. (2.10) versus Eq. (2.11)).

In either case, we denote by $\bar{A}_L$ the apparent conductivity of $A_L$ on $Q_L$ endowed with periodic boundary conditions, i.e., in analogy to the stochastic setting (2.1), for any $\bar{\xi} \in R^d$, via spatial averaging

$$ \bar{A}_L \bar{\xi} = \int_{Q_L} A_L \xi_L dx, \qquad (2.4) $$

instead of ensemble averaging and where $\xi_L \in L^2(Q_L; R^d)$ solves the equations

$$ \xi_L = \bar{\bar{\xi}} + \Gamma_L \xi_L \quad \text{and} \quad \Gamma_L A_L \xi_L = 0 \qquad (2.5) $$

in terms of the Helmholtz projector $\Gamma_L$ onto gradient fields of periodic functions on the cube $Q_L$. Notice that the apparent conductivity tensor $\bar{A}_L$ is still a random variable. Also, its expectation $\langle \bar{A}_L \rangle$ will be different from the effective conductivity $\bar{A}$,
<figure><img src="image_6.png" /><figcaption>Fig. 3. Periodizing the distance for a two-dimensional cell $Q_L$ containing two circles (red and green) of equal radius and their shaded periodically replicated images.</figcaption></figure>

in general, which is well-known in the engineering community, see Huet (1990) and Sab (1992) for early accounts. However, it can be shown, either in case $A_L = A_L^{sn}$ (Sab, 1992; Bourgeat and Piatnitski, 2004; Owhadi, 2003) or in case $A_L = A_L^{per}$, see Gloria et al. (2015) and Clozeau et al. (2021) for specific frameworks, that

$$ \bar{A}_L \to \bar{A} \quad \text{as} \quad L \to \infty \qquad (2.6) $$

with probability one, i.e., the apparent conductivity tensors $\bar{A}_L$ approximate the effective conductivity tensor $\bar{A}$ for increasingly large cubes $Q_L$. In the engineering community, a sufficiently large cell $Q_L$ is called a representative volume element (Hill, 1963; Kanit et al., 2003). Note that boundary conditions different from periodic could be used as well for defining the truncated corrector problem (2.5), e.g., Dirichlet or Neumann boundary conditions, but we chose periodic boundary conditions as they tend to be more accurate for fixed cell size, see, for instance, Kanit et al. (2003).

The convergence statement (2.6) is purely qualitative and does not permit to analyze the effect of the size of the volume elements. Indeed, in practice a suitable trade-off between the increased accuracy that accompanies larger computational cells and the limited computational resources needs to be identified, for instance based on statistical estimation techniques (Kanit et al., 2003). For this purpose, quantitative convergence statements become useful. On average, the error may be decomposed in the form

$$ \langle \|\bar{A}_L - \bar{A}\|^2 \rangle = \langle \|\bar{A}_L - \langle \bar{A}_L \rangle \|^2 \rangle + \langle \|\bar{A}_L\| - \bar{A}\|^2 \rangle \qquad (2.7) $$

in terms of the Frobenius norm, i.e., the right-hand side consists of a random error, which measures the standard deviation of the snapshot ensemble, and a systematic error which quantifies the defect induced by working on a cell of finite size. From a physical point of view, the random error quantifies the lack of statistical representativity of the realization in the finite-sized cell, whereas the systematic error accounts for artificial long-range correlations induced by periodization for $A_L^{per}$ and the boundary-layer error in case of taking snapshots $A_L^{sn}$.

By using a Monte-Carlo sampling-strategy on cells of fixed size, the random contribution may be decreased — in contrast to the systematic error. Indeed, suppose that $N$ independent coefficient fields ($A_{L,i}$) are sampled. Then, the empirical expectation of the corresponding effective conductivities satisfies

$$ \langle \left\| \frac{1}{N} \sum_{i=1}^{N} \bar{A}_{L,i} - \bar{A} \right\|^2 \rangle = \frac{1}{N} \left\langle \left\| \bar{A}_L - \langle \bar{A}_L \rangle \right\|^2 \right\rangle + \left\langle \langle \bar{A}_L \rangle - \bar{A} \right\|^2. $$

The decomposition (2.7) is implicit in the statistical approach to estimating the minimum size of a representative volume element pioneered by Kanit et al. (2003). They assume the systematic error, which they call bias, to be negligible compared to the random error, which they refer to as dispersion.

In general, due to the central limit theorem (CLT), we cannot expect a decay of the random error better than

$$ \sqrt{\langle \|\bar{A}_L - \langle \bar{A}_L \rangle \|^2 \rangle} \lesssim L^{-\frac{d}{2}}. \qquad (2.8) $$

This scaling was confirmed by Kanit et al. (2003) in computational experiments for synthetic Poisson–Voronoi microstructures and thermal conductivity. In case of elasticity, a slightly inferior scaling was observed.
On the theoretical side, notice that the decay of the error (2.7) may be arbitrarily slow for general stationary and just qualitatively ergodic coefficient fields. Specific quantified ergodicity assumptions permit drawing stronger conclusions. For instance, ensembles with finite range of dependence were considered by Armstrong and Smart (2016), who provided quantitative stochastic homogenization results for convex integral functionals of quadratic growth. Armstrong and Mourrat (2016) relaxed the finite-range assumption and provided results for a mixing condition with algebraic decorrelation rate, previously considered by Yurinskii (1986).

As an alternative to the approaches just mentioned, ergodicity may be quantified in terms of suitable functional inequalities in probability (Naddaf and Spencer, 1998), which compensate the lack of “natural” Poincaré inequality on the probability space. Under an $L^2$ spectral gap assumption, either if $A_L$ is obtained as a snapshot (Egloffe et al., 2015, eq. (3.3)) or by periodization (Duerinckx et al., 2020, Th. 1) (Gloria et al., 2015, Th. 2), the CLT estimate (2.8) is expected to hold. For the latter case, Duerinckx et al. (Duerinckx et al., 2020, Th. 1) actually identify the first order term in the expansion of $\bar{A}_L^{\text{per}} - \langle \bar{A}_L^{\text{per}} \rangle$. As a consequence, the CLT scaling (2.8) of the random error is indeed optimal.

However, these quantifications of ergodicity only hold under restrictive assumptions on the underlying random coefficient field and may not meet the requirements of practical models of interest to the applied sciences (Torquato, 2002). Thus, spectral gap assumptions allowing for thicker stochastic tails were considered, i.e., in terms of an $L^p$ spectral gap assumption (Fischer and Otto, 2017) and via a weighted logarithmic Sobolev inequality (Gloria et al., 2019, Def. 1). Further weighted functional inequalities satisfied by ensembles of relevance to materials science were studied by Duerinckx and Gloria (2020b,a). Note that the CLT scaling (2.8) may not be valid for all materials, as in the case of Poisson fibers with infinite length (Jeulin, 2016; Dirrenberger et al., 2014), which feature an inferior scaling. Also, in specific situations, a faster decay rate than the CLT scaling is possible. As a specific example in one spatial dimension, let us consider an ensemble with microstructures made up of rods with length $L_0$, where one half is occupied by a conducting medium with conductivity $a_1$ and the other half features a conductivity $a_2$. Then, each microstructure realization is made up of a series of connected, possibly flipped rods, where the flipping probability is one half. The resulting ensemble is stationary and ergodic, and the effective (scalar) conductivity $\bar{A}$ is just the harmonic mean of $a_1$ and $a_2$. Moreover, the apparent conductivity of a random sample of length $L$ has a conductivity $\bar{A}_L$ which satisfies

$$ \frac{|\bar{A} - \bar{A}_L|}{\bar{A}} \le \frac{L_0}{2L} \max(a_1, a_2) \left| \frac{1}{a_1} - \frac{1}{a_2} \right| \quad \text{for } L \ge L_0, \qquad (2.9) $$

which is faster than the CLT scaling $L^{-\frac{1}{2}}$ in this dimension. This medium, whose heat flux enjoys a vanishing integral range (Lan-tuéjoul, 1991, p. 393), is, however, non-generic. Concerning the systematic error, the convergence rate is sensitive to the way the finite-cell ensemble $A_L$ is built. In general, a boundary-layer error (as introduced earlier) with dimension-independent scaling

$$ \| \langle \bar{A}_L^{\text{sn}} \rangle - \bar{A} \| \lesssim L^{-1} \qquad (2.10) $$

is expected, see Egloffe et al. (2015, Eq. (3.4)), and proved to be optimal in Clozeau et al. (2021) for a specific class of coefficient fields. If the coefficient fields $A_L^{\text{sn}}$ are considered as given, techniques based on modifying the corrector Eq. (2.5) were proposed (Mourrat, 2019; Abdulle et al., 2019) to screen the influence of the improperly treated boundary of the considered cell and to restore a better convergence rate. In contrast, using the periodized coefficient field $A_L = A_L^{\text{per}}$ may attenuate the systematic error, so that screening techniques are not necessary.

The latter statement could be made rigorous for specific classes of ensembles that allow for a natural periodization. For i.i.d. discretely random conductivities, Gloria et al. (Gloria et al., 2015, Th. 2) established the estimate

$$ \| \langle \bar{A}_L^{\text{per}} \rangle - \bar{A} \| \lesssim L^{-d} \qquad (2.11) $$

for $d > 2$, where there is an additional factor of $\ln^d(L)$ on the r.h.s. Such a result is known to hold for the case of overlapping Poisson random inclusions. For the case where the ensemble is generated from a Gaussian field with integrable covariance function, Clozeau et al. (2021) showed estimate (2.11), and proved this scaling to be optimal in this setting by identifying the first order term of the asymptotic expansion of $\langle \bar{A}_L^{\text{per}} \rangle - \bar{A}$.

### 3. Computational tools

In this Section we provide details on the microstructure-generation methods used as the basis for the computational experiments conducted in Section 4.

#### 3.1. The mechanical contraction method

For composites with spherical fillers, the simple random sequential adsorption algorithm of Feder (1980) is unable to reach industrial-scale volume fractions in reasonable time. Indeed, the jamming limit of mono-sized spheres is approximately 38% (Meakin, 1992), which is insufficient for the computational experiments of Section 4.1 below, taking into account the isolation distance.

For this purpose, collective rearrangement algorithms are necessary, and several methods are available (Lubachevsky and Stillinger, 1990b; Torquato and Jiao, 2010). We shall describe the mechanical-contraction method (MCM) of Williams and Philipse (2003), as there is a modification which generalizes to cylindrical fillers. In its original form, the MCM was designed for generating dense packings, but we use it for generating packings of spheres with lower filler fraction, as well. The MCM consists of two different phases which are executed in an alternated fashion: contraction and overlap removal. Upon contraction, a configuration of non-overlapping spheres in a cell $Q_L$ is modified by shrinking the entire cell by a specific factor $\rho \in (0, 1)$, keeping the sphere's radii
<figure><img src="image_7.png" /><figcaption>Fig. 4. Illustration of the overlap-removal technique central to the mechanical contraction method (Williams and Philipse, 2003), applied to 25 circles and an area fraction of 80%.</figcaption></figure>

unchanged, but modifying the sphere centers via a $\rho$-scaling w.r.t. the background Cartesian coordinate system. If the contraction procedure leads to a configuration without overlap, the reached volume fraction will be increased by a factor of $1/\rho^3$. In case of overlap, the second phase of the MCM comes into play — the overlap removal. Suppose that $N$ spheres with radius $r$ and centers $x_1, \dots, x_N$ are given in the cubic cell $Q_L$. The MCM applies an overlap-removal technique in order to move the centers $x_1, \dots, x_N$ in such a way that the spheres do not overlap anymore. For this purpose, an overlap energy

$$W : Q_L^N \to \mathbb{R}, \quad (x_1, \dots, x_N) \mapsto \frac{1}{2} \sum_{1 \le i < j \le N} \delta(x_i, x_j)^2, \qquad (3.1)$$

is defined in terms of an overlap indicator

$$\delta(x_i, x_j) = \langle 2r - \mathrm{dist}(x_i, x_j) \rangle_+, \qquad (3.2)$$

where $\mathrm{dist}(x_i, x_j)$ denotes the $Q_L$-periodic distance of $x_i$ and $x_j$ and

$$\langle z \rangle_+ = \max(0, z)$$

is the Macaulay bracket.

Two (open) spheres with radius $r$ centered at $x_i$ and $x_j$, respectively, do not overlap if and only if $\delta(x_i, x_j) = 0$. Thus, the $N$ spheres centered at $x_1, \dots, x_N$ are in a non-overlapping configuration precisely if $W(x_1, \dots, x_N) = 0$. As the overlap energy $W$ (3.1) is continuously differentiable, a gradient-descent method may be used for finding a global minimizer of $W$, see Williams and Philipse (2003) for details. In dimension two, the approach is illustrated in Fig. 4, where 25 disks were placed in an initially random configuration at 80% area fraction. During the iterations, the centers of overlapping disks are moved, whereas disks without overlap remain fixed. For this example, the overlap-removal algorithm converged to the desired precision after 288 iterations. The MCM algorithm used in this work operates as follows. Suppose a decreasing sequence $L_1, L_2, \dots, L_M$ of volume-element edge-lengths is given, together with a fixed radius $r$ and an initial configuration $x_1, x_2, \dots, x_N$ of sphere centers, contained in $Q_{L_1}$. If the initial configuration has overlaps, the overlap-removal technique is applied. Then, for any $m = 1, \dots, M-1$, the sphere centers are shrunk by a factor $L_{k+1}/L_k$, followed by the overlap-removal stage. If the final overlap removal was successful, a sphere packing with a volume fraction of $\phi = N \times \frac{4}{3}\pi r^3/L_M^3$ is reached.
Some care has to be taken, as there are limits on the volume fractions to be reached. Also, the number of substeps M as well as the size of the shrinking factors $L_{k+1}/L_k$ have an influence on the achievable filler fractions. For the study at hand, we fix both the radius r and the volume fraction $\phi$, and suppose that the number of spheres N is prescribed. Then, we use the prescription

$$L_m = \frac{m}{M} r \sqrt[3]{\frac{4\pi N}{3\phi}}, \quad m = 1, 2, \dots, M,$$

with M = 3.

Up to this point, we described the deterministic algorithm where the number of spheres and the initial configurations are fixed. To create the ensembles, whose study is of primary concern to this article, this data may be drawn from suitable probability distributions. Typically, we consider the number of spheres to be fixed. Only in Section 5, we consider the case where the sphere count N is Poisson distributed. Once the number of spheres is fixed, we sample the initial centers of the spheres from a uniform distribution on the cube $Q_{L_1}$. We close describing the MCM algorithm with a few remarks.

1. The described strategy works for dimension two with obvious modifications.

2. The computational effort of the MCM algorithm is hidden in computing the pairwise distances (3.2). A naive implementation requires on the order of $N^2$ computations, which can be excessive for large N. We rely upon cell-linked lists (Mattson and Rice, 1999) to reduce the complexity to $O(N)$.

3. In practice, an isolation distance between the spheres is useful to avoid singularities in the solution field. This is realized by working with a slightly larger radius $\bar{r} > r$ when computing the overlap indicator (3.2).

## 3.2. Sequential addition and migration

The mechanical contraction method of Williams and Philipse (2003) also applies to spherocylinders, i.e., cylinders with spherical caps attached, provided the overlap indicator (3.2) is modified. A microstructure generated in this way may be used as a model of a fiber-reinforced composite by neglecting the spherical caps of the spherocylinders.

In its original form, however, the method can only generate moderate volume fractions without altering the overall fiber orientation. For this reason, Schneider (2017) modified the overlap energy W (3.1) by a penalty term accounting for the orientation. Also, in order to reach higher volume fractions, instead of successively shrinking the cell, fibers are added incrementally to a cell of fixed size. The resulting algorithm is called sequential addition and migration (SAM), see Schneider (2017) for details about the deterministic algorithm.

To specify the studied ensembles, it is necessary to quantify the uncertainty of the initial data of the SAM algorithm. In this study, we fix the number N of cylindrical fibers. The centers of the cylinders are drawn from a uniform distribution on the considered cube, whereas the long axes of the cylinders follow a uniform distribution on the unit sphere. The latter is drawn by sampling a vector from a three-dimensional standard normal distribution which is subsequently normalized.

# 4. Computational results

## 4.1. Setup

We are concerned with microstructures of particle-reinforced composites. We used the isotropic thermal conductivity of polypropylene (0.2 W/(m K)) for the matrix and of E-glass 1.2 W/(m K)) for the fillers, see Weidenfeller et al. (2005).

The microstructures were discretized on a regular pixel/voxel grid, segmented according to which phase the center of the pixel/voxel belongs to. Regular voxel meshes (Kim and Swan, 2003; Lian et al., 2013; Talebi et al., 2019) are commonplace for computational homogenization of industrial-scale microstructures, despite recent advances in mesh generation (Zhao et al., 2016; Sohn, 2018; Ullah et al., 2019). The static thermal equilibrium problem (2.5) was discretized by the Moulinec-Suquet discretization (Moulinec and Suquet, 1994, 1998), an underintegrated Fourier-Galerkin discretization (Vondřejc et al., 2014), see Brisard and Dormieux (2012) and Schneider (2015) for convergence proofs. For resolving the discrete system, we used the Eyre-Milton solver (Eyre and Milton, 1999), a Douglas-Rachford type solution method (Schneider et al., 2019), with a tolerance of $10^{-6}$, implemented in Python with Cython extensions, see Schneider (2019) for more details. We use periodic boundary conditions for the snapshots, although the microstructures are not periodic. Periodic boundary conditions come naturally with the Moulinec-Suquet discretization (Moulinec and Suquet, 1994, 1998).

For the circular and spherical inclusions, a desktop computer with a 6-core Intel i7 CPU and 32GB RAM was utilized. A dual-socket workstation with 48 cores and 1 TB RAM was used for the short-fiber microstructures.

## 4.2. Spherical inclusions

We investigate mono-disperse spherical fillers as our first example. For generating microstructures, we follow two protocols. The first strategy, which we call periodized, generates $N = K^3$ ($K = 2, 4, 8, 16$) non-overlapping spheres by the mechanical contraction method, see Section 3.1. For the second, snapshot strategy, $N = (2K)^3$ non-overlapping spheres are packed as before, but we cut out 1/8th of the structure. Thus, we arrive at a non-periodic structure of identical size as for the previous protocol. For the first protocol, $N = K^3$ spheres are chosen to minimize artifacts for small computational cells. Indeed, a regular configuration would support $K^3$
<figure><img src="image_8.png" /><figcaption>Fig. 5. 3D comparison of periodically generated cells and increasing number of inclusions, with $2^3 = 8$, $4^3 = 64$, $8^3 = 512$ and $16^3 = 4096$ spheres (from left to right).</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 6. A large volume element with $K = 64$, i.e., $64^3 = 262144$ spheres used for computing the reference effective conductivity (4.2).</figcaption></figure>

spheres, whereas integers close to $K^3$ may lead to artificial gaps in the structure and increase the systematic error. Please notice also that the first protocol does not involve any stochastic fluctuation when it comes to the achieved volume fraction. Indeed, by fixing both the filler count and the volume element, the reached volume fraction is fixed. In practice, however, our voxel-based meshing strategy introduces a (small) fluctuation in the volume fraction. In contrast, the second protocol will cut out $N = K^3$ spheres on average and a variation of the volume fraction is present. Slightly abusing terminology, we parameterize the cell size by the number of fillers $N = K^3$, also for the second protocol.

For both approaches, we fill to 30% by the mechanical contraction method of Williams and Philipse (2003) in 10% steps, and we use an isolation distance of 20% of the sphere's radius. Results of this microstructure-generation process are shown in Fig. 5.

A unit cell containing $K^3$ spheres (on average) is discretized by $(16K)^3$ voxels, i.e., the microstructures for $K = 16$ are discretized by $256^3 = 16777216$ elements and contain 4096 spheres. Following Khoromskia et al. (2020), Section 5, we generated 10000 microstructures for $K = 2, 4, 8, 16$ and the two protocols. The most time-consuming step here was generating the snapshot elements with $K = 16$, because 10000 volume elements, each containing 32768 non-overlapping spheres, needed to be generated quickly.
<caption>**Table 1**</caption>

Computed effective conductivities $\bar{a}_L$ (left, with 99% confidence intervals) and standard deviations (right) in W/(m K) for the spheres, computed using 10000 realizations and ten realizations for $K = 64$.

<table><thead><tr><th>K</th><th>Periodized</th><th>Snapshots</th><th>K</th><th>Periodized</th><th>Snapshots</th></tr></thead><tbody><tr><td>2</td><td>0.345553 ± 5 × 10<sup>-5</sup></td><td>0.356449 ± 4 × 10<sup>-4</sup></td><td>2</td><td>0.0017786</td><td>0.0172593</td></tr><tr><td>4</td><td>0.345538 ± 2 × 10<sup>-5</sup></td><td>0.351227 ± 2 × 10<sup>-4</sup></td><td>4</td><td>0.0006614</td><td>0.0080232</td></tr><tr><td>8</td><td>0.345382 ± 6 × 10<sup>-6</sup></td><td>0.348215 ± 1 × 10<sup>-4</sup></td><td>8</td><td>0.0002450</td><td>0.0043529</td></tr><tr><td>16</td><td>0.345291 ± 2 × 10<sup>-6</sup></td><td>0.346698 ± 5 × 10<sup>-5</sup></td><td>16</td><td>0.0000871</td><td>0.0020374</td></tr><tr><td>64</td><td>0.345224 ± 1 × 10<sup>-5</sup></td><td>–</td><td>64</td><td>0.0000107</td><td>–</td></tr></tbody></table>

<figure><img src="image_10.png" /><figcaption>**Fig. 7.** Convergence behavior of the systematic and the random error for spherical inclusions and the periodized/snapshot protocols, normalized according to Eq. (4.2).</figcaption></figure>

<caption>**Table 2**</caption>

Empirical probability (in %, 10000 realizations) of being 1%-close (left) and 0.1%-close (right) to $\bar{a}$ (4.2) for spherical fillers and 30% volume fraction.

<table><thead><tr><th>K</th><th>Periodized</th><th>Snapshots</th><th>K</th><th>Periodized</th><th>Snapshots</th></tr></thead><tbody><tr><td>2</td><td>95.77</td><td>12.20</td><td>2</td><td>17.36</td><td>1.42</td></tr><tr><td>4</td><td>100.00</td><td>24.06</td><td>4</td><td>35.29</td><td>2.35</td></tr><tr><td>8</td><td>100.00</td><td>45.45</td><td>8</td><td>76.06</td><td>4.68</td></tr><tr><td>16</td><td>100.00</td><td>82.73</td><td>16</td><td>99.93</td><td>10.13</td></tr></tbody></table>

For each of the 10000 generated microstructures, we computed the 11-component of the effective thermal conductivity. To keep notation simple, we introduce

$$ \bar{a} = \bar{A}_{11} \qquad (4.1) $$

as our quantity of interest. In contrast to the volume fraction, the ground truth for $\bar{a}$ is not known exactly. For this reason, we need to resort to computations on a sufficiently large scale. More precisely, we generated ten large microstructures with the periodization strategy, see Fig. 6, with $64^3 = 262144$ spherical inclusions, discretized by $1024^3 \approx 1.07 \times 10^9$ voxels, and computed the reference

$$ \bar{a} = 0.345228 \pm 0.000015 \text{ W/(m K)}. \qquad (4.2) $$

Here, Student's t-distribution for $t = 10$ is used for estimating a 99% two-sided confidence interval based on the (-1)-shifted standard deviation computed for these ten samples. Compared to the lower Hashin-Shtrikman bound (Hashin and Shtrikman, 1962), $\bar{a}^{HS} \approx 0.338462 \text{ W/(m K)}$, the computed value $\bar{a}$ is about 2% higher. Table 1 contains the computed means and standard deviations for $K = 2, 4, 8, 16$ and the periodized/snapshot protocols. The mean values are supplemented by 99%-two sided confidence intervals, computed using the standard deviations, and rounded to the highest significant digit. Due to the large number of samples used, these confidence intervals are much tighter for $8 \le K \le 16$ than for the ten computations on large volume elements, see Fig. 6.

The systematic error for the periodized ensembles is about an order of magnitude smaller than for their snapshot counterparts, see Fig. 7(a). For the latter, an $L^{-1}$-convergence behavior is evident. From an engineering point of view, a systematic error below 0.1% is sufficient — and this is computed using only eight (!) inclusions.

Taking a look at the relative standard deviation, see Fig. 7(b), we observe that the snapshot-sampling strategy leads to a convergence rate of $L^{-1}$. In contrast, the standard deviation for the periodized ensemble decreases as $L^{-3/2}$. Also, there is a difference by an order of magnitude between the two protocols.
<figure><img src="image_11.png" /><figcaption>Fig. 8. First of the ten microstructures with 1024² disks used for computing the effective conductivity (4.3).</figcaption></figure>

Due to proper scaling of the systematic and random errors in Fig. 7, we may also quantify the uncertainty involved in the computation. Indeed, the random error is about an order of magnitude larger than the systematic error for both protocols.

Mean and standard deviations only provide a limited amount of information concerning the full distribution of the random variable $a_L$. For this purpose, and from an engineering perspective, we may ask the following question: Suppose we run only a single computation, what is the probability of being 1% (0.1%) close to the reference (4.2)?

Using the 10000 computations for each setup, the empirical probabilities are collected in Table 2. We see that, for the periodized protocol, more than two correct digits are computed with a chance of more than 95% using $2^3 = 8$ spheres. For higher $K$, each of the 10000 computations predicted two significant digits correctly. In contrast, for the snapshot protocol, even using $16^3 = 4096$ spheres led to a lower success probability than for eight spheres in the periodized setting.

Computing the third significant digit with a single run is more challenging, and $K \ge 8$ is required for the periodized protocol to have a higher success than failure probability. Only for $K = 16$, i.e., 4096 spheres, three significant digits may be computed almost surely. Indeed, only 7 out of 10000 runs failed to provide three significant digits.

In contrast, for the snapshot approach, the success probability is about an order of magnitude smaller.

### 4.3. Circular inclusions

Subsequent to studying spherical inclusions, we turn to circular inclusions in two spatial dimensions. Microstructures of this type serve as models for continuous fibers, i.e., parallel cylindrical inclusions whose length is much larger than the diameter, and the effective behavior may be determined from a two-dimensional model. In turn, we may investigate the changes in convergence rates for stochastic homogenization when decreasing the ambient dimension.

We followed a similar protocol as for the spherical case, with the notable difference that we only need to cut out 1/4th from a sample containing $(2K)^2$ disks to obtain the snapshot. To ensure some degree of compatibility, we retained the area-fraction steps and the isolation distance for the MCM algorithm. We also used $(16K)^2$ pixels for discretizing the geometries filled by $K^2$ disks, on average. Due to the two-dimensional situation, we were able to treat larger $K$ than for the three-dimensional case.

For reference, we generated ten microstructures with a periodized protocol with $1024^2 = 1048576$ circular inclusions, see Fig. 8, discretized by $16384^2$ pixels. We obtained a reference conductivity, measured transverse to the fibers, of

$$ \bar{a} = 0.3174257 \pm 0.0000214 \text{ W/(m K)}, \qquad (4.3) $$

including a two-sided 99%-confidence interval based on Student's t-distribution for the ten drawn samples, and using identical material parameters as for the spherical inclusions. Notice that the effective conductivity (4.3) transverse to continuous fibers is smaller than the effective conductivity of spherical fillers (4.2) at the same filler fraction. Compared to the (two-dimensional) lower Hashin–Shtrikman bound (Hashin and Shtrikman, 1962), $\bar{a}^{HS} \approx 0.309090 \text{ W/(m K)}$ the computed value $\bar{a}$ is about 2.6% higher.

To assess the convergence behavior of the transverse effective conductivity, we generated 10000 microstructures for $K$ ranging from 2 to 64 in dyadic steps, both for the periodized and the snapshot protocol, and computed the $\bar{a}$-value. The results are collected in Table 3, together with the estimated standard deviations. For comparison, also the computed value for $K = 1024$ was included. As for the spherical case, the confidence interval for the $K = 1024$ reference computation is comparatively large,
<table><caption>**Table 3**</caption><caption>Computed effective conductivities $\bar{a}_L$ (left, with 99% confidence intervals) and standard deviations (right) in W/(m K) for the circular inclusions, computed using 10000 realizations and ten realizations for K = 1024.</caption><thead><tr><th>K</th><th>Periodized</th><th>Snapshots</th><th>K</th><th>Periodized</th><th>Snapshots</th></tr></thead><tbody><tr><td>2</td><td>0.316286 ± 2 × 10<sup>-4</sup></td><td>0.324902 ± 1 × 10<sup>-3</sup></td><td>2</td><td>0.007254</td><td>0.045318</td></tr><tr><td>4</td><td>0.317640 ± 1 × 10<sup>-4</sup></td><td>0.320024 ± 7 × 10<sup>-4</sup></td><td>4</td><td>0.004890</td><td>0.026951</td></tr><tr><td>8</td><td>0.317517 ± 6 × 10<sup>-5</sup></td><td>0.318881 ± 4 × 10<sup>-4</sup></td><td>8</td><td>0.002419</td><td>0.014757</td></tr><tr><td>16</td><td>0.317505 ± 3 × 10<sup>-5</sup></td><td>0.318172 ± 2 × 10<sup>-4</sup></td><td>16</td><td>0.001181</td><td>0.007808</td></tr><tr><td>32</td><td>0.317459 ± 2 × 10<sup>-5</sup></td><td>0.317617 ± 1 × 10<sup>-4</sup></td><td>32</td><td>0.000592</td><td>0.004039</td></tr><tr><td>64</td><td>0.317438 ± 7 × 10<sup>-6</sup></td><td>0.317591 ± 5 × 10<sup>-5</sup></td><td>64</td><td>0.000291</td><td>0.002059</td></tr><tr><td>1024</td><td>0.317426 ± 2 × 10<sup>-5</sup></td><td>—</td><td>1024</td><td>0.000021</td><td>—</td></tr></tbody></table>

<figure><img src="image_12.png" /><figcaption>**Fig. 9.** Convergence behavior of the systematic and the random error for circular inclusions and the periodized/snapshot protocols, normalized according to Eq. (4.3).</figcaption></figure>

<table><caption>**Table 4**</caption><caption>Empirical probability (in %, 10000 realizations) of being 1%-close (left) and 0.1%-close (right) to $\bar{a}$ (4.3) for circular inclusions.</caption><thead><tr><th>K</th><th>Periodized</th><th>Snapshots</th><th>K</th><th>Periodized</th><th>Snapshots</th></tr></thead><tbody><tr><td>2</td><td>32.35</td><td>5.61</td><td>2</td><td>3.47</td><td>0.67</td></tr><tr><td>4</td><td>50.27</td><td>9.47</td><td>4</td><td>5.47</td><td>0.92</td></tr><tr><td>8</td><td>80.93</td><td>17.52</td><td>8</td><td>10.81</td><td>1.78</td></tr><tr><td>16</td><td>99.25</td><td>30.95</td><td>16</td><td>21.30</td><td>3.24</td></tr><tr><td>32</td><td>100.00</td><td>56.77</td><td>32</td><td>40.60</td><td>5.86</td></tr><tr><td>64</td><td>100.00</td><td>87.47</td><td>64</td><td>72.71</td><td>12.31</td></tr></tbody></table>

approximately identical to the K = 32 case. This effect is a consequence of using 10 samples for K = 1024 instead of 10000 samples for K = 32. Also notice that K = 64 corresponds to K² = 4096 circular inclusions, discretized on a regular grid with 1024² pixels.

We observe, see Fig. 9(a), that the systematic error of the periodized protocol is consistently smaller than for the snapshot approach, almost by an order of magnitude. For the snapshot ensemble, the systematic error follows an L⁻¹-rate, as predicted. Evaluating the periodized protocol is much more difficult, as the systematic error is extremely small. In particular, and in view of the confidence interval of the reference computation, see Table 3, we cannot safely provide a convergence rate. For the standard deviation, see Fig. 9(b), both protocols lead to an L⁻¹-convergence rate. Thus, in contrast to the spherical case, there is no difference in convergence rates between the periodized and the snapshot ensemble. We will discuss the reasons behind this effect in Section 5. The standard deviation for the snapshot approach is about an order of magnitude larger than for the periodized setting. Comparing the levels of stochastic and systematic error, see Fig. 9, for both approaches the random error is about an order of magnitude larger than the systematic error.

As for the spherical inclusions treated in Section 4.2, we take a look at the empirical success probabilities of obtaining two and three correct significant digits, respectively, with only a single computation, see Table 4. For the periodized protocol, 16² circular inclusions are required to obtain two correct digits with 99% probability, i.e., more than 100 fibers are required. In contrast, for the snapshot strategy, 32² = 1024 inclusions are required to exceed 50% success probability, a value that is exceeded for 4² = 16 inclusions in case of the periodized protocol.

Obtaining three correct significant digits is difficult in the two-dimensional setting, even for the periodized protocol. Indeed, even working with 64² = 4096 inclusions leads to a failure rate of about 30%.
<figure><img src="image_13.png" /><figcaption>Fig. 10. Periodically generated short-fiber reinforced elements with increasing fiber count.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 11. Large volume element with $L/ℓ = 10$, including 76397 fibers. Smaller cell sizes with $L/ℓ \in \{1, 2, 4\}$ are indicated by white boxes in the lower left corner.</figcaption></figure>

### 4.4. Short-fiber reinforced composites

Last but not least, we consider short-fiber reinforced composites, i.e., microstructures consisting of non-overlapping cylindrical inclusions. More precisely, we consider identically shaped cylindrical fibers with length $ℓ$ and diameter $d$ at an aspect ratio $ℓ/d$ of 20, with an isotropic fiber orientation, encoded by an isotropic fourth-order fiber-orientation tensor (Advani and Tucker, 1987). In the periodized setting, we fill to 15% volume in three 5%-steps by Sequential Addition and Migration (Schneider, 2017), see Section 3.2, with a minimum distance of 20% of the fiber's diameter, and consider cubic volume elements $Q_L$, s.t. the edge length $L$ is an integer multiple of the fiber length $ℓ$, i.e., $L/ℓ = 1, 2, 3, 4$, see Fig. 10 for examples. We refer to Altendorf et al. (2014) for a related study with the periodized protocol.

For the snapshot protocol, we generate cells $Q_{3/2L}$ with the periodized protocol and extract a subcell of dimension $L^3$. As in the previous two Sections, we endow the matrix with the (isotropic) conductivity parameters of PP (0.2 W/(m K)) and the fibers with those of E-glass (1.2 W/(m K)). Each fiber is discretized by 5 voxels per diameter (Müller et al., 2015), i.e., 100 voxels per length $ℓ$. This translates into $(L/ℓ \times 100)^3$ voxels in total, i.e., $100^3$ voxels for $L = ℓ$ and $400^3$ voxels for $L = 4ℓ$.
<caption>**Table 5**</caption>

Computed effective conductivities $\bar{a}_L$ (left, with 99% confidence intervals) and standard deviations (right) in W/(m K) for short fibers, computed using 10000 realizations and ten realizations for $L/\ell = 10$.

<table><thead><tr><th>L/ℓ</th><th>Periodized</th><th>Snapshots</th><th>L/ℓ</th><th>Periodized</th><th>Snapshots</th></tr></thead><tbody><tr><td>1</td><td>0.281079 ± 6 × 10<sup>-6</sup></td><td>0.280321 ± 1 × 10<sup>-4</sup></td><td>1</td><td>0.000223</td><td>0.003678</td></tr><tr><td>2</td><td>0.281211 ± 2 × 10<sup>-6</sup></td><td>0.280583 ± 5 × 10<sup>-5</sup></td><td>2</td><td>0.000087</td><td>0.001823</td></tr><tr><td>3</td><td>0.281060 ± 1 × 10<sup>-6</sup></td><td>0.280716 ± 3 × 10<sup>-5</sup></td><td>3</td><td>0.000047</td><td>0.001178</td></tr><tr><td>4</td><td>0.281042 ± 1 × 10<sup>-6</sup></td><td>0.280766 ± 2 × 10<sup>-5</sup></td><td>4</td><td>0.000031</td><td>0.000837</td></tr><tr><td>10</td><td>0.281008 ± 6 × 10<sup>-6</sup></td><td>–</td><td>10</td><td>0.000006</td><td>–</td></tr></tbody></table>

<figure><img src="image_15.png" /><figcaption>**Fig. 12.** Convergence behavior of the systematic and the random error for short fibers and the periodized/snapshot protocols, normalized according to Eq. (4.4).</figcaption></figure>

<caption>**Table 6**</caption>

Empirical probability (in %, 10000 realizations) of being 1%-close (left) and 0.1%-close (right) to $\bar{a}$ (4.2) for fibers.

<table><thead><tr><th>L/ℓ</th><th>Periodized</th><th>Snapshots</th><th>L/ℓ</th><th>Periodized</th><th>Snapshots</th></tr></thead><tbody><tr><td>1</td><td>100.00</td><td>55.09</td><td>1</td><td>77.12</td><td>6.14</td></tr><tr><td>2</td><td>100.00</td><td>86.54</td><td>2</td><td>81.91</td><td>11.63</td></tr><tr><td>3</td><td>100.00</td><td>98.02</td><td>3</td><td>100.00</td><td>17.85</td></tr><tr><td>4</td><td>100.00</td><td>99.91</td><td>4</td><td>100.00</td><td>24.59</td></tr></tbody></table>

We use the average of 10 computed effective conductivities for samples with $L/\ell = 10$ (1000³ voxels, see Fig. 11) as the reference:

$$ \bar{a} = 0.281008 \pm 0.0000055 \text{ W/(m K)}, \qquad (4.4) $$

including a 99% two-sided confidence interval, as before. Compared to the lower Hashin–Shtrikman bound (Hashin and Shtrikman, 1962), $\bar{a}^{HS} \approx 0.262069 \text{ W/(m K)}$, the computed value $\bar{a}$ is about 6.7% higher. The empirical means of 10000 computed effective conductivities for the considered volume-element sizes and the two protocols are listed in Table 5, together with two-sided 99% confidence intervals and standard deviations. For completeness, the reference data (4.4) is also included.

Notice that the confidence intervals are very tight in this setting. Indeed, due to the shape of the inclusions, even for $L = \ell$, twenty fiber diameters fit across an edge of such a volume element. In fact, the confidence interval for $L = \ell$ is very close to the confidence interval for $L = 10\ell$, as 10000 samples were evaluated for the former, and only ten samples for the latter. Even for the snapshot-sampling strategy the random errors are small.

This accuracy is also reflected in the normalized systematic and random errors, see Fig. 12(b). For all settings considered, they do not exceed 2%.

For the systematic error, see Fig. 12(a), the snapshot strategy shows the already familiar $L^{-1}$-scaling, whereas the periodized sampling protocol appears to scale as $L^{-3}$ for $L \ge 2\ell$. Between $L = \ell$ and $L = 3\ell$, the systematic error does not decrease monotonically, presumably due to a geometric influence. Indeed, for $L = \ell$, an edge of the volume element under consideration is exactly equal to a fiber length. The comparatively low systematic error for this case may thus be a results of the severe restrictions in fiber arrangement for $L = \ell$. In magnitude, the systematic error for the snapshot approach is about half an order of magnitude larger than for the periodized sampling.

For the random error, see Fig. 12(b), we observe an $L^{-1}$ and $L^{-3/2}$ scaling for the snapshot and periodized protocol, respectively, in agreement with the observations made for spheres, see Fig. 7. Also, there is about an order of magnitude difference in the absolute
<figure><img src="image_16.png" /><figcaption>Fig. 13. Convergence behavior of the standard deviation for computing the volume fraction $\phi_L$ and periodized/snapshot protocols.</figcaption></figure>

standard deviations for both approaches, already evident in Table 5. It is interesting to note that, for the periodized sampling, the random error is on the same order of magnitude as the systematic error. This behavior contrasts what we observed for spherical and circular inclusions.

As for the other inclusion types, the empirical probabilities for remaining within 1% (and 0.1%) relative error for only a single computation are given in Table 6.

For the periodized sampling strategy, each of the 40000 computation was correct to two significant digits. Also, the chance for obtaining three correct digits exceeds 3/4 for $L = \ell$, already. For $L = 3\ell$ and above, each individual computation was guaranteed to feature three correct digits. In contrast, for the snapshot sampling, the success probabilities are much smaller. An edge length $L = 3\ell$ is required to obtain two correct digits with 98% probability, and the involved effort is at least $3^3 = 27$ times as high as using an element from the periodized ensemble and $L = \ell$. Obtaining higher accuracy than two digits appears out of reach for the snapshot strategy. Even for $L = 4\ell$, only a quarter of the 10000 computations led to results with three correct leading digits.

## 5. On the scaling of the random error in three dimensions

This section is intended to complement the results of the previous section. More precisely, we provide insights into the difference in the decay behavior of the random error for the snapshot ensembles compared to their periodized counterparts. It will also become clear why we observe this behavior in three spatial dimensions only.

For a random conducting medium with two distinct (positive) conductivities $\alpha_1$ and $\alpha_2$, a series expansion in terms of the material-contrast related quantity

$$ \rho = \frac{\sqrt{\alpha_1} - \sqrt{\alpha_2}}{\sqrt{\alpha_1} + \sqrt{\alpha_2}} $$

permits identifying the leading-order term in the estimate

$$ \sqrt{\langle \| \bar{A} - \bar{A}_L \|^2 \rangle} \le 4 \sqrt{\alpha_1 \alpha_2} \sqrt{\mathrm{Var}(\phi_L)} \rho + O(\rho^2). \quad (5.1) $$

We refer to the Appendix for a self-contained derivation. Thus, the decay rate of the random error is closely tied to the decay rate of the variance of the sampled volume fraction. In Fig. 13(a), we compare the standard deviations of the volume fractions $\phi_L$ for the periodized and the snapshot protocol and spherical fillers.

As the number of inclusions is fixed to $N = K^3$ for the periodized ensemble, the volume fractions are exact up to voxelation, i.e., the rastering process which furnishes each voxel either by matrix or filler material. The latter process, as it is based on the voxel's centroid, shows an $L^{-3/2}$-scaling in accordance with the CLT. For the snapshot ensemble, however, we observe an $L^{-1}$-scaling. In addition, we also sampled the number of spherical inclusions from a Poisson distribution, and used the MCM to generate a corresponding “periodized Poisson ensemble”. Unfortunately, if $K$ was too small, the number of cells produced by the Poisson process often led to a filler content that could not be realized by MCM. For this reason, only results for $K \ge 4$ are shown. Fig. 13(a) shows an $L^{-3/2}$-decay of the standard deviation for the periodized Poisson ensemble. However, the absolute values are on a higher level than for the previously discussed protocols. Note that the mean volume fractions for all three considered scenarios and all considered $K$ were within 0.1% relative error of the target 30% volume fraction.
<figure><img src="image_17.png" /><figcaption>(a) Single realizations at $L/\ell = 32$ and varying fiber (b) Mean and standard deviation for ten realizations at volume content $\phi$ $L/\ell = 64$ and $\phi = 30\%$</figcaption></figure>

<caption>Fig. 14. Magnitude of the autocorrelation function $|h|$, depending on $\|x\|$.</caption>

For comparison, the decay rates of the sampled volume fractions of the circular inclusions are shown in Fig. 13(b). Both the periodized and the snapshot-sampling strategies are characterized by a $L^{-1}$-scaling, in accordance with the CLT.

As an intermediate result, we conclude that the non-standard scaling of the random error, which we observed in Section 4, is already present at the level of the volume fraction. Fortunately, understanding the sampling of the volume fraction is simpler than understanding the random error through homogenization. Indeed, there is a (more) direct link between the autocorrelation function and the decay behavior of the volume fraction for the snapshot protocol, which arises as the mean of the (random) characteristic-function snapshot on cells $Q_L$ of finite size.

For estimating the decay rate of the volume fraction for snapshot ensembles, the notion of integral range (Matheron, 1989; Lantuëjoul, 2002) is very useful, and we will quickly summarize the important points.

Suppose that $Z$ is a stationary random field on $\mathbb{R}^d$ with finite second moment. By stationarity, the mean

$$ \mu = \langle Z(x) \rangle \quad \text{and the variance} \quad \sigma^2 = \mathrm{Var}(Z(x)) \equiv \langle (Z(x) - \mu)^2 \rangle $$

are independent of $x \in \mathbb{R}^d$. For any compact set $V \subset \mathbb{R}^d$, let us denote by $Z_V$ the random variable obtained by averaging $Z$ over $V$

$$ Z_V = \int_V Z(x) dx. $$

$Z_V$ gives rise to an unbiased estimator for the mean $\mu$. For $\sigma > 0$, we are interested in the variance of $Z_V$,

$$ \mathrm{Var}(Z_V) = \left\langle \left( \int_V Z(x) dx - \mu \right)^2 \right\rangle \equiv \sigma^2 \int_V \int_V h(x-y) dx dy, \qquad (5.2) $$

expressed in terms of the scaled autocorrelation function

$$ h(x) = \sigma^{-2} \langle (Z(0) - \mu)(Z(x) - \mu) \rangle. \qquad (5.3) $$

The identity (5.2) uncovers the relation between the variance of the $V$-averaged variable $Z_V$ and the autocorrelation function (5.3). If $h \in L^1$, the estimate

$$ \mathrm{Var}(Z_V) \le \sigma^2 \int_V \int_V |h(x-y)| dx dy \le \frac{\sigma^2}{|V|} \int_V \int_{\mathbb{R}^d} |h(x-y)| dx dy = \frac{\sigma^2 \|h\|_{L^1}}{|V|}, $$

holds, i.e., when evaluated on the cube $V = Q_L$,

$$ \mathrm{Var}(Z_{Q_L}) \le \sigma^2 \|h\|_{L^1} L^{-d}. $$

In particular, $Z_V$ converges to the mean $\mu$ (at least) with the CLT scaling. Thus, if the autocorrelation function $h$ decays to zero sufficiently rapidly, a CLT scaling is ensured. For slower decorrelations, inferior convergence rates may be expected.

For the problem at hand, we are interested in $Z = \chi$, the characteristic function of the random inclusions. Then, we calculate

$$ \mu = \phi \quad \text{and} \quad \sigma^2 = \phi(1 - \phi) $$

in terms of the particle volume fraction $\phi$. Also, in terms of the identification $\phi_L^{sn} = \chi_{Q_L}$, Eq. (5.2) becomes

$$ \mathrm{Var}(\phi_L^{sn}) = \phi(1 - \phi) \int_{Q_L} \int_{Q_L} h(x - y) dx dy. \qquad (5.4) $$
<figure><img src="image_18.png" /><figcaption>Fig. 15. Normalized standard deviation of the volume fraction and the effective conductivity for spherical particles.</figcaption></figure>

For spherical inclusions, we investigate the autocorrelation function h more closely, see Fig. 14.

Indeed, if the autocorrelation function h has a sufficiently rapid decay, the volume fraction $\phi_L$ of the snapshot will decay (at least) with the CLT scaling. Conversely, a non-CLT scaling must be rooted in slowly decaying decorrelations.

We considered cells with $32^3$ spheres at $\phi = 30\%$ as our point of departure, fixed their size and reduced the generated volume fractions. The resulting empirical autocorrelation functions, for a single realization and computed by an FFT-based approach (Fullwood et al., 2008), are shown in Fig. 14(a). The autocorrelation function h at x provides the probability of finding an inclusion both at y and y + x (averaged over y). As we consider only statistically isotropic ensembles, the autocorrelation function depends only on the magnitude $||x||$ of the vector x. Thus, the oscillations in Fig. 14(a) indicate closely packed spheres with a minimum distance between them. This is a result of the enforced minimum distance between the inclusions. With this interpretation in mind, Fig. 14(a) may be inspected more closely. For increasing volume fraction, the number of “links” in closely packed “chains” of spheres increases. Indeed, even for $\phi = 20\%$, only two spheres will be found in the vicinity of a fixed sphere. The third next sphere will be further apart already. For $\phi = 30\%$, already seven spheres are found close to a given sphere. For comparison, we also included the $\phi = 35\%$ case, which is even more highly packed. Actually, it was not possible to generate $\phi = 40\%$ with the MCM method. Indeed, due to the 20% isolation distance, we generate packings of spheres with a larger radius, and down-scale those spheres only in post-processing. For 20% isolation distance and a volume fraction $\phi = 40\%$, we would be generating sphere packings without isolation distance and a volume fraction $\phi = 1.2^3 \times 40\% = 69.12\%$. This number exceeds the jamming limit for spheres in three dimensions, and is thus extremely difficult, if not possible, to be realized by a random packing.

To assess the reliability of the computed autocorrelation functions, we computed ten realizations of the empirical autocorrelation functions for $64^3$ spheres at 30% volume fraction, and considered the mean and standard deviation of these samples, see Fig. 14(b). Notice that the mean autocorrelation function reaches the level of the standard deviation for about ten sphere radii and at a level of about $10^{-3}$. Thus, the actual decay rate of the autocorrelation function cannot be identified from these numerical investigations.

Thus, we turn our attention to Eq. (5.4), and assess the normalized quantity

$$ \sqrt{\frac{\mathrm{Var}(\phi_L^{\mathrm{sn}})}{\phi(1-\phi)}} = \sqrt{\int_{Q_L} \int_{Q_L} h(x-y) \, dx \, dy.} \qquad (5.5) $$

directly.

For increasing volume fraction, the normalized standard deviation (5.5) is shown in Fig. 15(a), where we introduced a reference length scale $L_0$, s.t. at $\phi = 30\%$, exactly $(L/L_0)^d$ spheres fit inside a cell with volume $L^d$ to realize this volume fraction. Thus, $L/L_0$ is identical to the parameter K used previously, but also makes sense for other values of the filler fraction $\phi$. Note that we only ran 1000 simulations to cover the great variety of considered volume fractions.

Interestingly, up to a volume fraction of 20%, the normalized standard deviation is independent of the volume fraction and decays as predicted by the CLT. Thus, we observe a regime of rapid decorrelation. Increasing the volume fraction further, in particular from 25% to 30%, leads to an incremental change in the convergence rate, continuously changing from the CLT scaling to a $L^{-1}$-rate. In Section 4.2, we considered precisely this regime. Increasing the volume fraction even further retains the $L^{-1}$-rate, but leads to an overall decrease of the standard deviation. Thus, we observe a regime with slow decorrelation but decreased variance due to an increased degree of spatial order.

The corresponding empirically computed, normalized standard deviation of the effective conductivity is shown in Fig. 15(b). Similar to the volume fraction, different regimes emerge, corresponding to rapidly decorrelating ensembles which a favorable convergence rate (slightly inferior than CLT) and ensembles with long-range order, and an associated $L^{-1}$-decay.

Motivated by these observations, we turn to the two-dimensional case of circular inclusions.
<figure><img src="image_19.png" /><figcaption>Fig. 16. Normalized standard deviation of the area fraction for circular inclusions.</figcaption></figure>

For spherical inclusions, the observed decay rate deteriorates close to the maximum packing fraction (accounting for the minimum distance). Thus, we similarly increased the packing area fraction for the circular inclusions, see Fig. 16. Using the MCM algorithm, area fractions up to $\phi = 55\%$ could be reached. Thus, in contrast to the three-dimensional case, the area fraction $\phi = 30\%$ considered in Section 4.3 is not even close to the maximum packing fraction. Thus, the corresponding normalized standard deviation of the snapshot area-fractions decays with the CLT scaling. Only for much higher area fractions, at about $\phi = 50\%$, an inferior decay rate of $L^{-1/2}$ appears. Interestingly, in such cases, there is a transition from a $L^{-1/2}$-scaling up to $K = 16$ to the expected $L^{-1}$-scaling for larger cells.

It appears reasonable that such a transition to the standard CLT scaling should also appear for spherical and cylindrical fillers in three spatial dimension. However, the necessary size of the volume elements exceeds current computational capacities. Hence, we cannot observe it in this study.

These guesses on a CLT scaling for sufficiently large cells can also be backed up by theory. Indeed, Jeulin (2016) showed that codimension-k linear varieties randomly dispersed in $R^d$ may lead to an $L^{d-k}$-scaling of the variance. Willot (2017) demonstrated that, for cylindrical fibers at low volume fraction, two scaling regimes for the variance of the volume fraction are observed. Indeed, for small cells, the fibers appear infinitely long, and an $L^{-1}$-scaling arises. For sufficiently large cells, larger than about ten (Willot, 2017, Sec. 4.3) times the fiber length, the fibers are “small” compared to the volume, and the expected CLT scaling is recovered.

For the spherical and the circular fillers, a similar behavior is observed. Indeed, due to the microstructure generation process, the individual fillers interact in chains, or even networks of chains, see Fig. 4. Thus, one-dimensional structures arise. The average length of these structures depends on the target volume fraction, and, in turn, determines the decay behavior, see Fig. 15(a).

The situation is visualized in Fig. 17. We compare two microstructures with the same microstructure parameters, but generated by different algorithms. The vanilla RSA method, see Fig. 17(a), leads to a configuration with a much more homogeneous appearance. In contrast, the MCM tends to form clusters, see Fig. 17(b). This does not come unexpected. Recall that there is a well-defined jamming limit for RSA, which limits the expected volume fraction to be reached. At this jamming limit, the microstructures attain a rather homogeneous appearance, essentially due to a lack of non-occupied spots.

In contrast, the MCM is a *collective re-arrangement* algorithm which permits reaching much larger volume fractions than the RSA. It will also reach some sort of uniform jamming state, but at a much higher filler content.

Returning to the *real* material, of course, we cannot ab initio decide which microstructure generation method is suitable for the microstructured material at hand. However, for matrix-inclusion composites with industrial filler content, we expect collective re-arrangement algorithms to be more appropriate. Indeed, sequential algorithms like RSA place the inclusions in a successive fashion. In particular, the correlation length of the emerging ensemble is on the order of the inclusion size, independent of the filler fraction. At high filler content, and also motivated by the real production processes, the individual particles are expected to actually affect their neighbors. Indeed, for instance during injection molding, the filler particles interact with each other, both, in terms of hydrodynamic interactions, and in terms of direct particle-particle collisions. Thus, the correlation length of the produced microstructured materials actually depends on the filler content, and may cover a few particle sizes for high volume fraction.

To sum up, the observed inferior decay of the random error for the cut microstructures is a consequence of slowly decaying decorrelations between the particles in the ensemble. The latter are expected to emerge at high filler content, independently of the microstructure generation algorithm used. Still, homogenization theory applies as long as the considered volume elements are large compared to the correlation length. Furthermore, in three dimensions and for the considered ensembles, we were unable to work on sufficiently large volume elements to reach the regime of the classical CLT scaling.

## 6. Conclusion

This work was devoted to a comparison of two sampling strategies for representative volume elements of matrix-inclusion composites. More precisely, we investigated a sampling strategy based on periodizing the ensemble and a snapshot approach which
<figure><img src="image_20.png" /><figcaption>Fig. 17. Visual comparison of two cells with 64² = 4096 circular inclusions at φ = 30% area fraction, generated by two different algorithms.</figcaption></figure>

extracts the sample from a larger realization. The latter is typical for microstructures obtained via digital imaging techniques, as real heterogeneous materials are not periodic, in general. Periodized sampling, on the other hand, is quite natural (Sab and Nedjar, 2005) when generating microstructures of matrix-inclusion composites synthetically. Our curiosity was driven by recent mathematical results (Clozeau et al., 2021), which show – under very specific assumptions on the ensemble – that the periodization protocol may be more advantageous than the snapshot approach in terms of the accuracy of the computed apparent properties. From an engineering point of view, the benefits of periodization were already realized by Sab and Nedjar (2005), who drew attention to the disturbance of statistical homogeneity introduced by taking a snapshot.

In the engineering literature, it is typical to rely upon only one of the two approaches, or, at least, to compare the results for the generated microstructures to either computations on a large digital image or experimental data. The objective of the work at hand was to study both approaches carefully, and to compare the resulting effective properties in the spirit of the seminal work of Kanit et al. (2003).

To profit from the periodized ensemble, we utilized periodic boundary conditions for the computations, as well. Imposing Dirichlet or Neumann boundary conditions on rectangular samples of the periodized ensemble disturbs statistical homogeneity in a similar way as working with the snapshot ensemble. For the snapshot ensembles, we utilized periodic boundary conditions, as they are known to lie between the Dirichlet and Neumann results, and the associated apparent properties are typically closer to the effective properties (Kanit et al., 2003; Ostoja-Starzewski, 2006; Salmi et al., 2012).

We considered the simple setting of thermal conductivity and considered the material parameters of E-glass fillers embedded in a polypropylene matrix, investigating spherical inclusions as well as continuous and short isotropic cylindrical fibers. For the computational investigations, we relied upon a voxel discretization (Kim and Swan, 2003; Lian et al., 2013; Talebi et al., 2019). It is well-known that voxelation introduces an additional error compared to the continuous case (Müller et al., 2015; Lian et al., 2013). However, once the voxel mesh size is fixed, the discretized microstructures may be regarded, both for the periodized ensemble and the snapshot strategy, as realizations of a corresponding *discrete* ensemble, very much like the scenario considered by Gloria et al. (2015). In particular, the changes in the statistical homogeneity of the samples would be present for the discretized setting, as well. Another source of potential error is given by the threshold of the termination criterion (Schneider, 2021, Sec. 3.6) used for the computations. To assess its magnitude, we repeated the computational experiments in two spatial dimensions (Section 4.3) with a tolerance that was one order of magnitude below the previously considered value. Changes in the plots shown in Fig. 9 were not visible for the naked eye.

In the computational experiments, we considered 10000 realizations for all but the largest microstructures, as, in previous studies (Khoromskia et al., 2020; Khoromskia and Khoromskij, 2020), this number was reported as sufficient to have enough confidence in the obtained mean apparent properties to deduce precise convergence rates. For all three considered cases, both the systematic error and the random error for the periodized ensemble were much lower than their snapshot counterparts by roughly an order of magnitude. Furthermore, for the three-dimensional examples, we observed a different scaling of the random error for the snapshot strategy than for the periodization approach. In Section 5, we showed that this behavior does no come unexpected in view of pertinent results from the field of spatial statistics (Matheron, 1989, 1971), in particular related to the notion of integral range (Lantuéjoul, 2002). Indeed, at high filler fraction, the particles in such composites may actually be strongly correlated to their neighbors in terms of their spatial position, leading to comparatively long tails in the correlation, encoded in terms of an increased integral range. Below or close to the integral range, the inferior scaling of the random error emerges.

This work sheds light on the seeming discrepancy between RVE sizes encountered for particle-reinforced composites. Indeed, the snapshot-based approaches typically predict a larger RVE size than reported for generated microstructures. Indeed, it should be
clear that a larger variance is expected for snapshots inside a large periodic volume compared to the use of a periodic small volume. As an example, let us take a glimpse at two studies on short-fiber reinforced composites. Altendorf et al. (2014) conduct simulations on large periodic cells, and compute the variance of the volume fraction as well as the thermal and elastic fields in non-overlapping subdomains. This well-established strategy (Escoda et al., 2016; Prill et al., 2017; Gasnier et al., 2015; Abdallah et al., 2015; Willot et al., 2016) aims to minimize edge effects, and is expected to provide “better” results than a pure snapshot strategy. In contrast, Schneider (2017) works with periodized ensembles generated by the SAM algorithm. Both simulation studies report on isotropic short-fiber structures with an aspect ratio of about 30. The elastic properties of the fibers coincide, but Schneider considers a higher filler content (20% vs 15%) and a lower material contrast (matrix with 3.3 GPa vs 2 GPa). Then, Altendorf et al. (Altendorf et al., 2014, Tab. 3) report that about 3.3 fiber lengths are necessary for volume elements to produce elastic properties with less than 1% standard deviation, whereas Schneider (Schneider, 2017, Sec. 3.3.1) reports two fiber lengths to yield about 0.15% relative error compared to a larger-scale result. Although the scenarios are not fully comparable, it appears that small periodized volumes may be beneficial when solely effective properties are of interest and computational effort is to be minimized.

It would be interesting to understand whether the inferior accuracy associated to the snapshot strategy may be attenuated via appropriate countermeasures, for instance by statistical methods (Ohser and Mücklich, 2000) or modified correctors (Mourrat, 2019; Abdulle et al., 2019). Also, the investigations underlined the benefits of proper morphological modeling (Jeulin, 2012) and microstructure-generation tools (Bargmann et al., 2018), even if digital images are available. Indeed, relying upon periodized ensembles may lead to substantial savings in computational effort compared to the snapshot strategy. It should be noted, however, that it is not always apparent how to periodize a given ensemble. For instance, Jeulin and coworkers (Jeulin, 2016; Dirrenberger et al., 2014) and Sukiman et al. (2017) consider infinitely long Poisson fibers. A naive periodization of this ensemble on the torus would almost surely be space-filling. Still, for many particle-reinforced composites, producing periodized samples within a microstructure-generation framework is rather obvious, and the benefits may be exploited.

In view of the work of Kanit et al. (2003), we could confirm their observation that the systematic error for periodized ensembles does not exceed the random error, at least for the material contrast considered, so that the random error is a good indicator for the total RVE error (Jeulin, 2012, Sec. 5.5).

We chose to work with thermal conductivity and a small material contrast to show the high probabilities for failing to produce accurate predictions of effective quantities, emphasizing that already for this seemingly innocent scenario, using the naive snapshot strategy may be suboptimal. Investigating higher material contrast, effective permeability (Abdallah et al., 2015; Willot et al., 2016) and mechanical problems, in particular inelastic problems, is of immediate engineering interest. Groundwork was already laid by Kanit et al. (2003), who reported an inferior convergence rate for the random error associated to mechanical problems. Also, Kanit et al. (2003) argued that the random error scales as the material contrast, which remains to be confirmed by theoretical means.

From the viewpoint of mathematical analysis, this article is a natural continuation of previous work by Khoromskia et al. (2020), Khoromskia and Khoromskij (2020), where quantitative homogenization results for random media were confirmed by numerical simulations. Indeed, we go beyond the cited work by investigating microstructure models of higher complexity, as is required for engineering applications. In particular, we break into terrain previously uncharted from a theoretical point of view. The ensemble which we studied differs from the simpler one of independently drawing particles conditioned on no overlap. The former permits to reach higher volume fraction than the latter, in general, which is a desirable property for engineering applications. As the studied ensemble reveals a cross-over in the error scaling at a volume-fraction dependent intermediate length scale, it appears interesting to investigate this ensemble in more detail.

As quantitative stochastic homogenization crucially relies upon quantitative versions of ergodicity, typically in the form of functional inequalities in probability, it may be of interest to study which conditions hold for the ensembles considered in this work, see Duerinckx and Gloria (2020b,a) for investigations covering Poisson–Voronoi tessellations, as treated by Kanit et al. (2003), and random sequential adsorption. In a similar direction, it might be interesting to study the intermediate length scale below which the convergence rate of the random error is inferior for snapshots of matrix-inclusion composites. Modified corrector problems (Gloria and Habibi, 2016; Mourrat, 2019) aim to upgrade the decay of the systematic error to the level of the random error, but typically assume a CLT scaling of the random error. For the studied matrix-inclusion composites and the snapshot protocol, a strictly inferior decay of the random error was observed in an intermediate regime for high filler fraction and moderate cell sizes. It appears sensible to investigate whether those methods that modify the corrector equation (Gloria and Habibi, 2016; Abdulle et al., 2019) improve the decay of the random error for such ensembles, as well. As opposed to Khoromskia et al. (2020), we could not (yet) confirm the optimal convergence rate of the systematic error and properly periodized ensembles predicted by Clozeau et al. (2021). Computations on even larger cells could reveal the proper decay rate of the systematic error.

A cornerstone of the seminal paper of Kanit et al. (2003) is the concept of integral range in spatial statistics (Lantuéjoul, 2002; Matheron, 1989), which quantifies the constant in front of the CLT scaling of a second-order stationary random field. In their study of the homogenization commutator, Duerinckx et al. (2020) related the Q-tensor, introduced by Mourrat and coworkers (Gu and Mourrat, 2016; Mourrat and Otto, 2016), to the stochastic fluctuations of the apparent properties associated to RVEs. Due to the characterization (Duerinckx et al., 2020, Thm. 2) in terms of a limiting procedure, the Q-tensor may be regarded as an extension of the integral range to stochastic homogenization, and a closer investigation appears desirable.

CRediT authorship contribution statement

**Matti Schneider:** Software, Validation, Formal analysis, Investigation, Writing – original draft, Writing – review & editing,
Visualization. **Marc Josien:** Methodology, Formal analysis, Writing – review & editing. **Felix Otto:** Conceptualization, Formal
analysis, Writing – review & editing.
Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgments

M. Schneider was supported by the German Research Foundation (DFG) within the International Research Training Group “Integrated engineering of continuous–discontinuous long fiber reinforced polymer structures” (GRK 2078). M. Schneider thanks S. Gajek for computational support and H. Andrä for stimulating discussions. This work was initiated while M. Josien held a post-doc position at the Max Planck Institute for Mathematics in the Sciences (MPI MiS), Leipzig. We thank the anonymous reviewers for their constructive feedback.

Appendix. A series expansion for the effective thermal conductivity

In this appendix, we collect the arguments leading to the estimate (5.1)

$$
\sqrt{\langle \| \bar{A} - A_L \|^2 \rangle} \le 4 \sqrt{\alpha_1 \alpha_2} \sqrt{\mathrm{Var}(\phi_L)} \rho + O(\rho^2)
$$

for the convenience of the reader. Let us fix some positive definite reference thermal conductivity tensor A⁰. By algebraic manipulations, see Milton (Milton, 2002, Sec. 14.9), for fixed $\bar{\xi} \in R^d$, the random vector field $\xi$ solves the Eqs. (2.2) and (2.3) precisely if $p = (A + A^0)\xi$ solves the Eyre–Milton equation (Eyre and Milton, 1999)

$$
(\mathrm{Id} - Y Z^0) p = 2 A^0 \bar{\xi}, \quad (\mathrm{A}.1)
$$

involving the Helmholtz reflection $Y = \mathrm{Id} - 2\Gamma$ and the Cayley mapping $Z^0 = (A - A^0)(A + A^0)^{-1}$. As the reflection $Y$ is orthogonal and $Z^0$ is an $L^2$-contraction provided the coefficient field $A$ is essentially bounded and uniformly positive definite (Eyre and Milton, 1999), the polarization $p$ solving the Eyre–Milton Eq. (A.1) may be represented in terms of a Neumann series

$$
p = 2 \sum_{k=0}^{\infty} (YZ^0)^k A^0 \bar{\xi} \qquad (A.2)
$$

with corresponding expectation

$$
\langle p \rangle = 2 \sum_{k=0}^{\infty} \langle (YZ^0)^k A^0 \bar{\xi} \rangle .
$$

Taking into account $p = q + A^0 \xi$, we obtain

$$
\bar{A} \bar{\xi} + A^0 \bar{\xi} = 2 A^0 \bar{\xi} + 2 \langle Y Z^0 A^0 \bar{\xi} \rangle + 2 \sum_{k=2}^{\infty} \langle (Y Z^0)^k A^0 \bar{\xi} \rangle,
$$

i.e., using that the expectation vanishes on the image of the Helmholtz projector Γ,

$$
\bar{A}_{\bar{\xi}} = A^{0}_{\bar{\xi}} + 2 \langle Z^0 A^0 \bar{\xi} \rangle + 2 \sum_{k=2}^{\infty} \langle (YZ^0)^k A^0 \bar{\xi} \rangle, \quad (A.3)
$$

To proceed, we restrict our attention to two-phase isotropic composites, i.e., we suppose that the thermal conductivity tensor takes the form

$$
A = \alpha_1 \chi \mathrm{Id} + \alpha_2 (1 - \chi) \mathrm{Id} \tag{A.4}
$$

in terms of non-vanishing isotropic thermal conductivities $\alpha_1$ as well as $\alpha_2$ and the indicator function $\chi$ of a random set. Define $A^0 = \alpha_0 \mathrm{Id}$ with $\alpha_0 = \sqrt{\alpha_1 \alpha_2}$. Then,

$$
Z^0 = (A - A^0)(A + A^0)^{-1} = \frac{\alpha_1 - \alpha_0}{\alpha_1 + \alpha_0} \chi \mathrm{Id} + \frac{\alpha_2 - \alpha_0}{\alpha_2 + \alpha_0} (1 - \chi) \mathrm{Id} = \frac{\sqrt{\alpha_1} - \sqrt{\alpha_2}}{\sqrt{\alpha_1} + \sqrt{\alpha_2}} [\chi - (1 - \chi)] \mathrm{Id},
$$

i.e.,

$$
Z^0 = \rho(2\chi - 1)\mathrm{Id} \quad \text{with} \quad \rho = \frac{\sqrt{\alpha_1} - \sqrt{\alpha_2}}{\sqrt{\alpha_1} + \sqrt{\alpha_2}} \in (-1, 1).
$$

Inserting the latter expression into the Neumann series (A.3), we obtain

$$
\bar{A}_{\bar{\xi}} = \sqrt{\alpha_1 \alpha_2 \bar{\xi}} + 2 \sqrt{\alpha_1 \alpha_2} \rho (2 \langle \chi \rangle - 1) \bar{\xi} + 2 \sqrt{\alpha_1 \alpha_2} \sum_{k=2}^{\infty} \langle (YZ^0)^k \bar{\xi} \rangle,
$$
which we might also write in the form

$$ \bar{A}\bar{\xi} = \sqrt{\alpha_1\alpha_2(1 - 2\rho)\bar{\xi} + 4\sqrt{\alpha_1\alpha_2\rho}\langle\chi\rangle\bar{\xi} + O(\rho^2)}. $$

We recover the known result, see Milton (Milton, 2002, Ch. 14), that, to first order in $\rho$, the volume fraction $\phi \equiv \langle\chi\rangle$ determines the effective conductivity.

For the snapshot ensemble on the cell $Q_L$, an analogous argument yields the representation

$$ \bar{A}^{\text{sn}} = \sqrt{\alpha_1\alpha_2(1 - 2\rho)\bar{\xi} + 4\sqrt{\alpha_1\alpha_2\rho}\phi_L^{\text{sn}}\bar{\xi} + O(\rho^2)}, $$

where the constant in the Landau-O depends on $L$ and the snapshot volume-fraction $\phi_L^{\text{sn}}$ is a random variable determined by

$$ \phi_L^{\text{sn}} = \int_{Q_L} \chi \, dx. $$

The snapshot volume-fraction $\phi_L^{\text{sn}}$ leads to an unbiased estimator for the volume fraction $\phi$, i.e.,

$$ \langle \phi_L^{\text{sn}} \rangle = \phi $$

holds. In particular, the systematic error $\|\bar{A} - \bar{A}^{\text{sn}}\|$ is entirely determined by the $\rho^2$-term. In contrast, the formula for the difference

$$ \bar{A} - \bar{A}^{\text{sn}} = 4\sqrt{\alpha_1\alpha_2}(\phi - \phi_L^{\text{sn}})\rho \text{ Id} + O(\rho^2) $$

permits bounding the random error in the form

$$ \sqrt{\langle \| \bar{A} - \bar{A}^{\text{sn}} \|^2 \rangle} \le 4\sqrt{\alpha_1\alpha_2}\sqrt{\text{Var}(\phi_L^{\text{sn}})\rho + O(\rho^2)}. \quad (A.5) $$

In particular, the random error is, up to first order in the material contrast dependent parameter $\rho$, determined by the standard deviation of the snapshot volume-fraction $\phi_L^{\text{sn}}$. We close this Section with a remark: The expansion (5.1) also holds in case of the periodized ensemble $A_L^{\text{per}}$. If the latter is constructed in such a way that the associated volume fraction for specific $L$ is such that $\phi_L^{\text{per}} = \phi$ holds, the linear term in $\rho$ will vanish for the standard-deviation expansion (5.1), i.e., the random error may be expected to become smaller. The latter idea has been extended to higher-order terms by Le Bris et al. (2016) and analyzed by Fischer (2019).

## References

* Abdallah, B., Willot, F., Jeulin, D., 2015. Stokes flow through a Boolean model of spheres: Representative volume element. Transp. Porous Media 109 (3), 711-726.

* Abdallah, B., Willot, F., Jeulin, D., 2016. Morphological modelling of three-phase microstructures of anode layers using SEM images. J. Microsc. 263 (1), 51-63.

* Abdulle, A., Arjmand, D., Paganoni, E., 2019. Exponential decay of the resonance error in numerical homogenization via parabolic and elliptic cell problems. C. R. Math. 357 (6), 545-551.

* Advani, S.G., Tucker, C.L., 1987. The use of tensors to describe and predict fiber orientation in short fiber composites. J. Rheol. 31 (8), 751-784.

* Altendorf, H., Jeulin, D., Willot, F., 2014. Influence of the fiber geometry on the macroscopic elastic and thermal properties. Int. J. Solids Struct. 51 (23), 3807-3822.

* Arjmand, D., Runborg, O., 2016. A time dependent approach for removing the cell boundary error in elliptic homogenization problems. J. Comput. Phys. 341, 206-227.

* Armstrong, S.N., Mourrat, J.-C., 2016. Lipschitz Regularity for elliptic equations with random coefficients. Arch. Ration. Mech. Anal. 219, 255-346.

* Armstrong, S.N., Smart, C.K., 2016. Quantitative stochastic homogenization of convex integral functionals. Ann. Sci. de L'Ecole Norm. Super. 49 (2), 423-481.

* Babuska, I., 1973. Solution of interface problems by homogenization I. SIAM J. Math. Anal. 7, 603-634.

* Bargmann, S., Klusemann, B., Markmann, J., Schnabel, J.E., Schneider, K., Soyarslan, C., Wilmers, J., 2018. Generation of 3D representative volume elements for heterogeneous materials: A review. Prog. Mater. Sci. 96, 322-384.

* Blanc, X., Le Bris, C., 2010. Improving on computation of homogenized coefficients in the periodic and quasi-periodic settings. Netw. Heterog. Media 5 (1), 1-29.

* Bouchedjra, M., Boulemia, C., Amrouche, A., 2018. Determination of the RVE size for polycrystal metals to predict monotonic and cyclic elastoplastic behavior: statistical and numerical approach with new criteria. Eur. J. Mech. A Solids 72, 1-15.

* Bourgeat, A., Piatnitski, A., 2004. Approximations of effective coefficients in stochastic homogenization. Annales de l'Intitut H. Poincaré 40, 153-165.

* Brisard, S., Dormieux, L., 2012. Combining Galerkin approximation techniques with the principle of Hashin and Shtrikman to derive a new FFT-based numerical method for the homogenization of composites. Comput. Methods Appl. Mech. Engrg. 217-220, 197-212.

* Chen, L., Gu, B., Tao, J., Zhou, J., 2019. Modelling mesoporous alumina microstructure with 3D random models of platelets. Compos. Struct. 216, 279-289.

* Clozeau, N., Josien, M., Otto, F., Xu, Q., 2021. Bias in the representative volume element method: periodize the ensemble instead of its realizations. in preparation.

* De Giorgi, E., Spagnolo, S., 1973. Sulla convergenza degli integrali dell'energia per operatori ellittici del secondo ordine. Ann. Sc. Norm. Super. Pisa Cl. Sci. 8, 391-411.

* Dirrenberger, J., Forest, S., Jeulin, D., 2014. Towards gigantic RVE sizes for 3D stochastic fibrous networks. Int. J. Solids Struct. 51 (2), 359-376.

* Doškář, M., Zeman, J., Jarušková, D., Novák, J., 2018. Wang tiling aided statistical determination of the representative volume element size of random heterogeneous materials. Eur. J. Mech. A Solids 70, 280-295.

* Drugan, W.J., Willis, J.R., 1996. A micromechanics-based nonlocal constitutive equations and estimates of representative volume element size for elastic composites. J. Mech. Phys. Solids 44, 497-524.

* Duerinckx, M., Gloria, A., 2020a. Multiscale functional inequalities in probability: Concentration properties. ALEA Lat. Am. J. Probab. Math. Stat. 17, 133-157.

* Duerinckx, M., Gloria, A., 2020b. Weighted functional inequalities: Constructive approach. Ann. Henri Lebesgue 3, 825-872.

* Duerinckx, M., Gloria, A., Otto, F., 2020. The structure of fluctuations in stochastic homogenization. Comm. Math. Phys. 377, 259-306.

* Egloffe, A., Gloria, A., Mourrat, J.-C., Nguyen, T.N., 2015. Random walk in random environment, corrector equation and homogenized coefficients: from theory to numerics, back and forth. IMA J. Numer. Anal. 35 (2), 499-545.

* El Moumen, A., Kanit, T., Imad, A., 2021. Numerical evaluation of the representative volume element for random composites. Eur. J. Mech. A Solids 86, 104181.
El Moumen, A., Kanit, T., Imad, A., El Minor, H., 2015. Effect of reinforcement shape on physical properties and representative volume element of particles-reinforced composites: statistical and numerical approaches. *Mech. Mater.* 83, 1–16.

Escoda, J., Willot, F., Jeulin, D., Sanahuja, J., Toulemonde, C., 2016. Influence of the multiscale distribution of particles on elastic properties of concrete. *Internat. J. Enggr. Sci.* 98, 60–71.

Evesque, P., 2005. Fluctuations, correlation and representative elementary volume (REV) in granular materials. *Poudres Grains* 11, 6–17.

Eyre, D.J., Milton, G.W., 1999. A fast numerical scheme for computing the response of composites using grid refinement. *Eur. Phys. J. Appl. Phys.* 6 (1), 41–47.

Feder, J., 1980. Random sequential adsorption. *J. Theoret. Biol.* 87 (2), 237–254.

Figliuzzi, B., Jeulin, D., Faessel, M., Willot, F., Koishi, M., Kowatari, N., 2016. Modelling the microstructure and the viscoelastic behaviour of carbon black filled rubber materials from 3D simulations. *Tech. Mech.* 32 (1–2), 22–46.

Fischer, J., 2019. The choice of representative volumes in the approximation of effective properties of random materials. *Arch. Ration. Mech. Anal.* 234 (2), 635–726.

Fischer, J., Otto, F., 2017. Sublinear growth of the corrector in stochastic homogenization: optimal stochastic estimates for slowly decaying correlations. *Stoch. Partial Differ. Equ. Anal. Comput.* 5, 220–255.

Fritzen, F., Forest, S., Böhlke, T., Kondo, D., Kanit, T., 2012. Computational homogenization of elasto-plastic porous metals. *Int. J. Plast.* 29, 102–119.

Fullwood, D.T., Niezgoda, S.R., Kalidindi, S.R., 2008. Microstructure reconstructions from 2-point statistics using phase-recovery algorithms. *Acta Mater.* 56 (5), 942–948.

Gasnier, J.B., Willot, F., Trumel, H., Figliuzzi, B., Jeulin, D., Biessy, M., 2015. A Fourier-based numerical homogenization tool for an explosive material. *Mater. Tech.* 103 (3), 308.

Ghossein, E., Lévesque, M., 2012. A fully automated numerical tool for a comprehensive validation of homogenization models and its application to spherical particles reinforced composites. *Int. J. Solids Struct.* 49 (11–12), 1387–1398.

Ghossein, E., Lévesque, M., 2015. Homogenization models for predicting local field statistics in ellipsoidal particles reinforced composites: Comparisons and validations. *Int. J. Solids Struct.* 58, 91–105.

Gitman, I.M., Askes, H., Sluys, L., 2007. Representative volume: Existence and size determination. *Eng. Fract. Mech.* 74 (16), 2518–2534.

Gitman, I.M., Gitman, M.B., Askes, H., 2006. Quantification of stochastically stable representative volumes for random heterogeneous materials. *Arch. Appl. Mech.* 75 (2–3), 79–92.

Gloria, A., 2008. An analytical framework for numerical homogenization. *II. Windowing and oversampling*. *Multiscale Model. Simul.* 7 (1), 274–293.

Gloria, A., 2011. Reduction of the resonance error - part 1: Approximation of homogenized coefficients. *Math. Models Methods Appl. Sci.* 21 (8), 1601–1630.

Gloria, A., Habibi, Z., 2016. Reduction in the resonance error in numerical homogenization II: Correctors and extrapolation. *Found. Comput. Math.* 16, 217–296.

Gloria, A., Neukamm, S., Otto, F., 2015. Quantification of ergodicity in stochastic homogenization: optimal bounds via spectral gap on glauber dynamics. *Invent. Math.* 199 (2), 455–515.

Gloria, A., Neukamm, S., Otto, F., 2019. Quantitative estimates in stochastic homogenization for correlated coefficient fields. pp. 1–39, *ArXiv:1910.05530*.

Gloria, A., Otto, F., 2016. Quantitative estimates on the periodic approximation of the corrector in stochastic homogenization. *ESAIM: Proc. Surv.* 48, 80–97.

Gu, Y., Mourrat, J.-C., 2016. Scaling limit of fluctuations in stochastic homogenization. *Multiscale Model. Simul.* 14 (1), 452–481.

Gusev, A.A., 1997. Representative volume element size for elastic composites: A numerical study. *J. Mech. Phys. Solids* 45 (9), 1449–1459.

Guven, I., Cinar, K., 2019. Micromechanical modeling of particulate-filled composites using micro-CT to create representative volume elements. *Int. J. Mech. Mater.* Des. 15, 695–714.

Harper, L.T., Qian, C., Turner, T.A., Li, S., Warrior, N.A., 2012. Representative volume elements for discontinuous carbon fibre composites. Part 2: Determining the critical size. *Compos. Sci. Technol.* 72 (2), 204–210.

Hashin, Z., Shtrikman, A., 1962. A variational approach to the theory of the effective magnetic permeability of multiphase material. *J. Appl. Phys.* 32, 3125–3131.

Hill, R., 1963. Elastic properties of reinforced solids: Some theoretical principles. *J. Mech. Phys. Solids* 11 (5), 357–372.

Hoang, T.H., Guerich, M., Yvonnet, J., 2016. Determining the size of RVE for nonlinear random composites in an incremental computational homogenization framework. *J. Eng. Mech.* 142 (5), 04016018.

Hou, T.Y., Wu, X.-H., 1997. A multiscale finite element method for elliptic problems in composite materials and porous media. *J. Comput. Phys.* 134, 169–189.

Huet, C., 1990. Application of variational concepts to size effects in elastic heterogeneous bodies. *J. Mech. Phys. Solids* 38, 813–841.

Jean, A., Jeulin, D., Forest, S., Cantournet, S., N'Guyen, F., 2011. A multiscale microstructure model of carbon black distribution in rubber. *J. Microsc.* 241 (3), 243–260.

Jeulin, D., 2012. Morphology and effective properties of multi-scale random sets: A review. *C. R. Mécanique* 340 (4), 219–229.

Jeulin, D., 2016. Power Laws Variance Scaling of Boolean Random Varieties. *Methodol. Comput. Appl. Probab.* 18 (4), 1065–1079.

Jiang, M., Ostoja-Starzewski, M., Jasiuk, I., 2001. Scale-dependent bounds on effective elastoplastic response of random composites. *J. Mech. Phys. Solids* 49 (3), 655–673.

Kanit, T., Forest, S., Galliet, I., Mounoury, V., Jeulin, D., 2003. Determination of the size of the representative volume element for random composites: statistical and numerical approach. *Int. J. Solids Struct.* 40 (13–14), 3647–3679.

Kanit, T., Forest, S., Jeulin, D., N'Guyen, F., Singleton, S., 2011. Virtual improvement of ice cream properties by computational homogenization of microstructures. *Mech. Res. Commun.* 38 (2), 136–140.

Kanit, T., N'Guyen, F., Forest, S., Jeulin, D., Reed, M., Singleton, S., 2006. Apparent and effective physical properties of heterogeneous materials: Representativity of samples of two materials from food industry. *Comput. Methods Appl. Mech. Engrg.* 195 (33–36), 3960–3982.

Khisaeva, Z., Ostoja-Starzewski, M., 2006. On the size of RVE in finite elasticity of random composites. *J. Elasticity* 85 (2), 153–173.

Khoromskaia, V., Khoromskij, B.N., 2020. Tensor-based techniques for fast discretization and solution of 3D elliptic equations with random coefficients. pp. 1–25, *ArXiv:2007.06524*.

Khoromskaia, V., Khoromskij, B.N., Otto, F., 2020. Numerical study in stochastic homogenization for elliptic partial differential equations: Convergence rate in the size of representative volume elements. *Numer. Linear Algebra Appl.* 27, e2296.

Kim, H.J., Swan, C.C., 2003. Voxel-based meshing and unit-cell analysis of textile composites. *Internat. J. Numer. Methods Engrg.* 56 (7), 977–1006.

Kouznetsova, V., Geers, M.G., Brekelmans, W.M., 2002. Multi-scale constitutive modelling of heterogeneous materials with a gradient-enhanced computational homogenization scheme. *Internat. J. Numer. Methods Engrg.* 54 (8), 1235–1260.

Kouznetsova, V., Geers, M.G., Brekelmans, W.M., 2004. Size of a representative volume element in a second-order computational homogenization framework. *Int. J. Multiscale Comput. Eng.* 2 (4), 575–598.

Kozlov, S.M., 1978. Averaging of differential operators with almost periodic rapidly oscillating coefficients. *Math. USSR Sbornik* 107 (149) (2 (10)), 199–217.

Landis, E.N., Keane, D.T., 2010. X-ray Microtomography. *Mater. Charact.* 61, 1305–1316.

Lantuéjoul, C., 1991. Ergodicity and integral range. *J. Microsc.* 161 (3), 387–403.

Lantuéjoul, C., 2002. Geostatistical Simulation - Models and Algorithms. Springer, Berlin.

Larsen, E., 1975. Neutron transport and diffusion in inhomogeneous media. *J. Math. Phys.* 16, 1421–1427.

Le Bris, C., Legoll, F., Minvielle, M., 2016. Special quasirandom structures: A selection approach for stochastic homogenization. *Monte Carlo Methods Appl.* 22 (1), 25–54.

Lian, W.D., Legrain, G., Cartraud, P., 2013. Image-based computational homogenization and localization: comparison between X-FEM/levelset and voxel-based approaches. *Comput. Mech.* 51, 279–293.
Lubachevsky, B.D., Stillinger, F.H., 1990b. Geometric properties of random disk packings. J. Stat. Phys. 60 (5-6), 561-583.

Ma, J., Sahraee, S., Wriggers, P., De Lorenzis, L., 2015. Stochastic multiscale homogenization analysis of heterogeneous materials under finite deformations with full uncertainty in the microstructure. Comput. Mech. 55 (5), 819-835.

Masson, D., Abdallah, B., Willot, F., Jeulin, D., Mercadelli, E., Sanson, A., Thorel, A., 2015. Morphological modelling of a metal foam supported SOFC configuration. ECS Trans. 68 (1), 2951-2960.

Matheron, G., 1971. The Theory of Regionalized Variables and Its Applications. Paris School of Mines Publication.

Matheron, G., 1989. Estimating and Choosing. Springer, Berlin.

Matouš, K., Geers, M.G.D., Kouznetsova, V.G., Gillman, A., 2017. A review of predictive nonlinear theories for multiscale modeling of heterogeneous materials. J. Comput. Phys. 330, 192-220.

Mattson, W., Rice, B.M., 1999. Near-neighbor calculations using a modified cell-linked list method. Comput. Phys. Comm. 119, 135-148.

Meakin, P., 1992. Random sequential adsorption of spheres of different sizes. Physica A 187 (3), 475-488.

Michel, J., Moulinec, H., Suquet, P., 1999. Effective properties of composite materials with periodic microstructure: a computational approach. Comput. Methods Appl. Mech. Engrg. 172 (1-4), 109-143.

Milton, G.W., 2002. The Theory of Composites. Cambridge University Press, Cambridge.

Ming, W.E.P., Zhang, P., 2005. Analysis of the heterogeneous multiscale method for elliptic homogenization problems. J. Amer. Math. Soc. 18 (1), 121-156.

Mirkhalaf, S.M., Pires, F.M. Andrade, Simoes, R., 2016. Determination of the size of the representative volume element (RVE) for the simulation of heterogeneous polymers at finite strains. Finite Elem. Anal. Des. 119, 30-44.

Moulinec, H., Suquet, P., 1994. A fast numerical method for computing the linear and nonlinear mechanical properties of composites. Comptes Rendus de L'Académie Des Sciences. Série II 318 (11), 1417-1423.

Moulinec, H., Suquet, P., 1998. A numerical method for computing the overall response of nonlinear composites with complex microstructure. Comput. Methods Appl. Mech. Engrg. 157, 69-94.

Mourrat, J., 2019. Efficient methods for the estimation of homogenized coefficients. Found. Comput. Math. 19, 435-483.

Mourrat, J.-C., Otto, F., 2016. Correlation structure of the corrector in stochastic homogenization. Ann. Probab. 44 (5), 3207-3233.

Moussaddy, H., Therriault, D., Lévesque, M., 2013. Assessment of existing and introduction of a new and robust efficient definition of the representative volume element. Int. J. Solids Struct. 50 (24), 3817-3828.

Müller, V., Kabel, M., Andrä, H., Böhlke, T., 2015. Homogenization of linear elastic properties of short fiber reinforced composites – a comparison of mean field and voxel-based methods. Int. J. Solids Struct. 67-68, 56-70.

Naddaf, A., Spencer, T., 1998. Estimates on the variance of some homogenization problems. pp. 1-14, Unpublished Preprint.

Nguyen, V.P., Lloberas-Valls, O., Stroeven, M., Sluys, L.J., 2010. On the existence of representative volumes for softening quasi-brittle materials-a failure zone averaging scheme. Comput. Methods Appl. Mech. Engrg. 199 (45), 3028-3038.

Niezgoda, S.R., Turner, D.M., Fullwood, D.T., Kalidindi, S.R., 2010. Optimized structure based representative volume element sets reflecting the ensemble-averaged 2-point statistics. Acta Mater. 58 (13), 4432-4445.

Ohser, J., Mücklich, F., 2000. Statistical Analysis of Microstructures in Materials Science. Wiley, New York.

Ostoja-Starzewski, M., 2006. Material spatial randomness: from statistical to representative volume element. Probab. Eng. Mech. 21 (2), 112-132.

Ostoja-Starzewski, M., Schulte, J., 1996. Bounding of effective thermal conductivities of multiscale materials by essential and natural boundary conditions. Phys. Rev. B 54 (1), 278.

Ostoja-Starzewski, M., Schulte, J., 1998. Random field models of heterogeneous materials. Int. J. Solids Struct. 35 (19), 2429-2455.

Owhadi, H., 2003. Approximation of the effective conductivity of ergodic media by periodization. Probab. Theory Related Fields 125, 225-258.

Papanicolaou, G.C., Varadhan, S.R.S., 1981. Boundary value problems with rapidly oscillating random coefficients. In: Random Fields, Vol. I, II (Esztergom, 1979). In: Colloq. Math. Soc. János Bolyai, vol. 27, North-Holland, Amsterdam-New York, pp. 835-873.

Pelissou, C., Baccou, J., Monerie, Y., Perales, F., 2009. Determination of the size of the representative volume element for random quasi-brittle composites. Int. J. Solids Struct. 46 (14-15), 2842-2855.

Prill, T., Jeulin, D., Willot, F., Balach, J., Soldera, F., 2017. Prediction of effective properties of porous carbon electrodes from a parametric 3D random morphological model. Transp. Porous Media 120 (1), 141-165.

Qidwai, S.M., Turner, D.M., Niezgoda, S.R., Lewis, A.C., Geltmacher, A.B., Rowenhorst, D.J., Kalidindi, S.R., 2012. Estimating the response of polycrystalline materials using sets of weighted statistical volume elements. Acta Mater. 60, 13-14.

Ranganathan, S.I., Ostoja-Starzewski, M., 2008. Scaling function, anisotropy and the size of RVE in elastic random polycrystals. J. Mech. Phys. Solids 56 (9), 2773-2791.

Sab, K., 1992. On the homogenization and the simulation of random materials. Eur. J. Mech. A Solids 11, 585-607.

Sab, K., Nedjar, B., 2005. Periodization of random media and representative volume element size for linear composites. C. R. Mécanique 333 (2), 187-195.

Salmi, M., Auslender, F., Bornert, M., Fogli, M., 2012. Various estimates of representative volume element sizes based on a statistical analysis of the apparent behavior of random linear composites. C. R. Mécanique 340 (4-5), 230-246.

Salmi, M., Auslender, F., Bornert, M., Fogli, M., 2012a. Apparent and effective mechanical properties of linear matrix-inclusion random composites: Improved bounds for the effective behavior. Int. J. Solids Struct. 49 (10), 1195-1211.

Saroukhani, S., Vafadari, R., Andersson, R., Larsson, F., Runesson, K., 2015. On statistical strain and stress energy bounds from homogenization and virtual testing. Eur. J. Mech. A Solids 51, 77-95.

Savvas, D., Stefanou, G., Papadrakakis, M., 2016. Determination of RVE size for random composites with local volume fraction variation. Comput. Methods Appl. Mech. Engrg. 305, 340-358.

Schneider, M., 2015. Convergence of FFT-based homogenization for strongly heterogeneous media. Math. Methods Appl. Sci. 38 (13), 2761-2778.

Schneider, M., 2017. The sequential addition and migration method to generate representative volume elements for the homogenization of short fiber reinforced plastics. Comput. Mech. 59, 247-263.

Schneider, M., 2019. On the Barzilai-Borwein basic scheme in FFT-based computational homogenization. Internat. J. Numer. Methods Engrg. 118 (8), 482-494.

Schneider, M., 2021. A review of non-linear FFT-based computational homogenization methods. Acta Mech. 232, 2051-2100.

Schneider, M., Wicht, D., Böhlke, T., 2019. On polarization-based schemes for the FFT-based computational homogenization of inelastic materials. Comput. Mech. 64 (4), 1073-1095.

Sebsadji, S.K., Chouicha, K., 2012. Determining periodic representative volumes of concrete mixtures based on the fractal analysis. Int. J. Solids Struct. 49 (21), 2941-2950.

Šedlbauer, D., Lepš, M., 2019. Wang tiling for particle heterogeneous materials: Algorithms for generation of tiles/cubes via molecular dynamics. Appl. Comput. Mech. 11 (1), 53-76.

Segurado, J., Llorca, J., 2002. A numerical approximation to the elastic properties of sphere-reinforced composites. J. Mech. Phys. Solids 50 (10), 2107-2121.

Shan, Z., Gokhale, A.M., 2002. Representative volume element for non-uniform micro-structure. Comput. Mater. Sci. 24 (3), 361-379.

Sohn, D., 2018. Periodic mesh generation and homogenization of inclusion-reinforced composites using an element-carving technique with local mesh refinement. Compos. Struct. 185, 65-80.

Stroeven, M., Askes, H., Sluys, L., 2004. Numerical determination of representative volumes for granular materials. Comput. Methods Appl. Mech. Engrg. 193 (30-32), 3221-3238.
Sukiman, M.S., Kanit, T., N'Guyen, F., Imag, A., El Moumen, A., Erchiqui, F., 2017. Effective thermal and mechanical properties of randomly oriented short and long fiber composites. *Mech. Mater.* 107, 56–70.

Swaminathan, S., Ghosh, S., 2006. Statistically equivalent representative volume elements for unidirectional composite microstructures: Part II-with interfacial debonding. *J. Compos. Mater.* 40 (7), 605–621.

Talebi, H., Silani, M., Klusemann, B., 2019. The scaled boundary finite element method for computational homogenization of heterogeneous media. *Internat. J. Numer. Methods Engrg.* 118 (1), 1–17.

Teferra, K., Graham-Brady, L., 2018. A random field-based method to estimate convergence of apparent properties in computational homogenization. *Comput. Methods Appl. Mech. Engrg.* 330, 253–270.

Terada, K., Hori, M., Kyoya, T., Kikuchi, N., 2000. Simulation of the multi-scale convergence in computational homogenization approaches. *Int. J. Solids Struct.* 37 (16), 2285–2311.

Torquato, S., 2002. Random Heterogeneous Materials – Microstructure and Macroscopic Properties. Springer, New York.

Torquato, S., Jiao, Y., 2010. Robust algorithm to generate a diverse class of dense disordered and ordered sphere packings via linear programming. *Phys. Rev. E* 82 (6), 061302.

Trias, D., Costa, J., Turon, A., Hurtado, J., 2006. Determination of the critical size of a statistical representative volume element (SRVE) for carbon reinforced polymers. *Acta Mater.* 54 (13), 3471–3484.

Ullah, Z., Zhou, X.-Y., Kaczmarczyk, L., Archer, E., McIlhagger, A., Harkin-Jones, E., 2019. A unified framework for the multi-scale computational homogenisation of 3D-textile composites. *Composites B* 167, 582–598.

Vondřejc, J., Zeman, J., Marek, I., 2014. An FFT-based Galerkin method for homogenization of periodic media. *Comput. Math. Appl.* 68 (3), 156–173.

Wang, H., Pietrasanta, A., Jeulin, D., Willot, F., Faessel, M., Sorbier, L., Moreaud, M., 2015. Modelling mesoporous alumina microstructure with 3D random models of platelets. *J. Microsc.* 260 (3), 287–301.

Weidenfeller, B., Höfer, M., Schilling, F.R., 2005. Thermal conductivity, thermal diffusivity, and specific heat capacity of particle filled polypropylene. *Composites A* 35 (4), 423–429.

Williams, S., Philipse, A., 2003. Random packings of spheres and spherocylinders simulated by mechanical contraction. *Phys. Rev. E* 67, 1–9.

Willot, F., 2017. Mean Covariogram of Cylinders and Applications to Boolean Random Sets. *J. Contemp. Math. Anal. (Arm. Acad. Sci.)* 52 (6), 305–315.

Willot, F., Abdallah, B., Jeulin, D., 2016. The Permeability of Boolean Sets of Cylinders. *Oil Gas Sci. Technol. Rev. DIFP Energies Nouvelles* 71 (4), 52.

Yue, X., W.E., 2007. The local microscale problem in the multiscale modeling of strongly heterogeneous media: effects of boundary conditions and cell size. *J. Comput. Phys.* 222 (2), 556–572.

Yurinskii, V.V., 1986. Averaging of symmetric diffusion in random medium. *Sib. Math. J.* 27, 603–613.

Zhao, J.-X., Coupez, T., Decencière, E., Jeulin, D., Cárdenas-Peña, D., Silva, L., 2016. Direct multiphase mesh generation from 3D images using anisotropic mesh adaptation and a redistancing equation. *Comput. Methods Appl. Mech. Engrg.* 309, 288–306.

Zhikov, V.V., Kozlov, S.M., Oleinik, O.A., 1994. Homogenization of Differential Operators and Integral Functionals. Springer-Verlag, New York.