Local elastic moduli of simple random composites computed at different length scales

Techniques like nanoindentation and atomic force microscopy can estimate the local elastic moduli in a region surrounding the probe used. For composites with phase regions much larger than the size of the probe, these procedures can identify the phases via their different elastic moduli but identifying phase regions that are on the same size scale as the indent is more problematic. This paper looks at three random 3D 8003 voxel composite models, each consisting of a matrix and spherical inclusions. One model has non-overlapping spheres and two models have overlapping spheres, with two and three distinct phases. The linear elastic problem is solved for each microstructure, and histograms are made of the local Young’s moduli over a number of sub-volumes (SVs), averaged over progressively larger SVs. The number and shape of histogram peaks change from N delta functions, where N is the number of elastically distinct phases, at the 1 voxel SV limit, to a single delta function located at the value of the effective global Young’s modulus, when the SV equals the unit cell volume. The phase volume fractions are also tracked for each bin in the Young’s modulus histograms, showing the phase make-up of bin in the histogram. There are clear differences seen between the non-overlapping and three-phase overlapping models and the two-phase overlapping sphere model, because of different size microstructural features, characterized by the average value of size as computed by the W(q) function. In the three-phase model, a peak that is originally all phase 3 persists at its same location, but as the size of the SVs increase, it is made up of a mixture of phases, so that it cannot be identified with a single phase even though it remains a clear peak. These results give some guidance as to what probe size might be useful in distinguishing different phases by local elastic moduli measurements, and how the length scales of the probe and the microstructure interact.


Introduction
For linear elastic random composite materials, the global effective moduli have long ago been defined [1][2][3][4][5]. One method to compute these moduli is to apply a strain tensor e o kl to a representative volume (RV) of composite material, solve the elastic equations, and compute the strain and stress everywhere in the volume, denoted e kl and r ij . An RV of a random composite material can be qualitatively defined as being a volume large enough so that the computed effective elastic moduli change only negligibly with increasing RV size. This size is a function primarily of the microstructure and secondarily of the elastic contrast between phases. The averages of these quantities over the complete volume are denoted \e kl [ and \r ij [. Equation (1) implicitly defines the effective elastic moduli tensor, C ijkl : Throughout this paper, however, the Voigt notation for the components of the stress, strain, and elastic moduli tensors is used: Equation (2) defines the effective moduli tensor. If periodic boundary conditions are used, as in this paper, the average strain tensor,\e j [, is equal to the applied strain tensor, e o j [1][2][3][4][5]. When comparing any kind of composite theory to experiment, the effective global moduli are typically used since most experiments measure the effective moduli, usually in some direction, through a simple stress/strain measurement. In the last few decades, techniques like nanoindentation [6] and atomic force microscopy (AFM) [7,8] can now estimate the local moduli in a region surrounding the probe. The size of these regions is thought to be a few cubic micrometers in volume for nanoindentation and significantly smaller for AFM. For composites with homogeneous phase regions much larger than the size of the probe, these procedures can identify different phase moduli and therefore estimate phase boundaries. For composites where the phase regions are on the same size scale or smaller than the probe size, identifying different phases is more difficult, since the elastic fields from the indenter can overlap two or more phases, mixing different phase properties into the elastic response of the probe. This ''mixing'' has been sometimes denoted the ''substrate effect'' when nanoindentation has been applied to composite materials [9]. This paper treats general composite material models but has in mind applications to composite materials like cement paste [8,9].
These developments lead to the problem of how to define the effective moduli at the local level, where Eq. (2) is defined for a sub-volume (SV) of the material, which is what is sensed by the kind of instruments mentioned above. How do these local moduli depend on the size of the SVs that define them? A problem similar in some ways to that studied in this paper was treated by Ostoja-Starzewski [10], who looked at computing elastic moduli over a sub-area of a 2D random composite (inclusions in a matrix) and studying how the effective moduli varied as a function of sub-area size. However, Ref. [10] only looked at single sub-areas of a given size, not multiple SVs dividing up an entire microstructure. Other workers have considered the related problem of the effect on computed elastic properties of the interplay between the size of the representative volume element and the internal length scale of microstructural features [11][12][13][14][15][16][17].
Imagine a material divided up into cubic voxels (see Fig. 1, left), where each voxel is a single phase, and a voxel is the smallest possible unit of volume; i.e. a digital image model. The elastic equations have been solved, and the average strain and stress tensor are known inside each voxel. If Eq. (2) is applied to compute the effective elastic moduli tensor in each voxel, there would be N distinct values, where N is the number of elastically distinct phases in the composite. If the volume fraction of material having a certain Young's modulus is plotted in histogram form, there would be N delta functions (in the infinite system limit) with the weight for each equal to the volume fraction of that particular phase. If, on the other hand, Eq. (2) were to be used over the entire volume (see Fig. 1, right), then the same sort of plot would be a delta function (in the infinite system limit), centered at the effective Young's modulus of the entire system with volume fraction weight equal to unity. The point of this paper is to see how those N delta functions, defined for single voxels, are transformed into a single delta function as the SVs over which Eq. (2) is defined increase from single voxels to the entire volume (see Fig. 1, middle). This transformation is a direct consequence of the competition between the characteristic length scale of individual phases and the length scale or size of the SV used.
With this general picture, one can then comment on the possibility of some kind of probe actually measuring true single-phase moduli when the measurement of effective moduli is taken over a certain length scale (volume). This paper looks at three simple 3D models, two overlapping sphere models and one non-overlapping sphere model, in a volume that is 800 3 voxels in size. The complete elastic problem is solved for each microstructure using a parallelized finite element method, and histograms are made of the local elastic moduli found by averaging over individual SVs. The phase volume fractions found in each SV are computed and averaged over all those SVs whose average E values contribute to a given bin in the elastic moduli histogram, so that one can tell directly how much of each phase contributes to various local elastic moduli peaks.

Finite element method and averaging procedure
A digital-image-based finite element method (FEM) designed to compute the elastic properties of random microstructures [18][19][20] was used in this paper. The code uses Message Passing Interface (MPI) parallelization [20] to enable computations on large systems, where each cubic voxel is a tri-linear finite element and the mesh is a 3D rectangular parallelepiped of cubic voxels. The nodes are just the corners of the cubic voxels. The left image in Fig. 1 can be thought of as a schematic picture of a typical 3D mesh for a cube of material. It has been successfully used to study simple shapes like the dilute limit of rectangular blocks [21,22] embedded in a periodic unit cell, with the aid of digital resolution scaling [23,24] to improve accuracy. The algorithm can also operate directly on the image stacks that come from techniques like X-ray computed tomography [25][26][27], after segmentation into distinct phases and where each image is thought of as a single layer of voxels.. The entire image stack used defines the 3D rectangular parallelepiped mesh. In the parallelization scheme, the system is divided into a number of multi-layer stacks oriented normal to the z direction. Each processor controls one stack that is several voxels thick and communicates with the processor below and above it. The typical number of processors used is nz/4 to nz/2, where nz = the total number of layers in the z direction, and nx and ny are the number of voxels in the x and y directions, respectively. The total number of degrees of freedom in this mesh is 3 nx ny nz, since there is one real node per voxel. Details of the parallelization process can be found in the manual for the parallel program [20]. The manual for the scalar version of the code gives more theoretical and computational details [19] of the finite element portion of the computation. Both the scalar and parallel versions are freely available on-line [28]. In this paper, the six different independent strains are applied, using periodic boundary conditions, and the resulting stresses and strains are averaged over each voxel and in principle stored. However, for 800 3 finite elements, this is a lot of storage, so there is a pre-averaging step, taken over every 4 3 voxel SV, 200 3 in all, and these averaged local stresses and strains are stored for each choice of applied strain. Not saving the stresses and Fig. 1 A schematic view of a representative volume element of a composite material with averaging schemes that start out subdivided into 1 voxel SVs (left), then larger SVs = collections of voxels (middle), and finally averaging the stress and strain over the entire material volume (right) strains for each voxel saves on storage requirements, and the results for 1 voxel SVs are analytically known, since applying the equivalent of Eq. (2) to the stresses and strains in a single voxel just recovers the assigned moduli for that voxel.
To define the local elastic moduli as a function of length scale means averaging the stress and strain over each SV of the entire microstructure and then determining the unique elastic moduli tensor for that SV. This is the equivalent of Eq. (2), but done for an SV, not over the entire 800 3 voxel/finite element unit cell. Since there are in general 36 independent components of the elastic moduli tensor, 36 linear equations are required to determine these components. These equations are obtained by using the stored average stresses and strains. Each strain component, e j , is applied one at a time and gives rise to six stress components connected to this strain component, r i , where i = 1,6 and r i is a linear function of the 36 C ij components. If this procedure is carried out six times, then 6 • 6 = 36 independent linear equations are generated for each SV, whose solution is the 36 components of the elastic moduli tensor defined for that SV. A simple Gaussian elimination algorithm is used to find the components of C ij for a given SV. The general symmetry elastic moduli tensor for each SV is then isotropically averaged [29] and the Young's modulus (E) extracted. All results in this paper are given for E, which is closest to the indentation modulus measured in nanoindentation; the bulk and shear moduli behaved very similarly and so these results are not shown. The average Poisson's ratio per SV can be similarly defined but will not be considered further here. This paper studies how the values of E, computed for each SV, are distributed statistically as a function of SV size. The phase volume fractions are computed in each SV, and then averaged over all the SVs whose E values contribute to a given bin in the E histogram. Therefore, the phase volume fractions are also shown in terms of histograms, showing how they are distributed as a function of SV size. Two limits are known: when the SV size goes to zero (one voxel, in this case) and when the SV size becomes equal to the entire unit cell size. In this latter limit, the average of the SVs must equal the average E of the entire system and the width of the single E peak narrows to a delta function, with weight 1, in the limit of a very large unit cell. In this delta function, the phase fractions take on their global values, since they have been averaged over the entire system. In the single voxel SV limit, there are N sharp peaks for N elastically distinct phases with weights being the phase fractions / i and located at E i . This paper looks at how these two distributions change into each other as the size of the SVs ranges between these two limits.

Microstructure models: generation and characterization
Two simple random microstructures were first constructed, with the cubic unit cell 800 voxels on a side [1], so that the total number of degrees of freedom in the finite element mesh will be 800 3 nodes times 3 DOF per node = 1.536 9 10 9 . This was large enough so that the volume fractions of the various phases in the various models to be discussed would not change appreciably between different realizations of the same microstructure. This allowed us to only use one realization of each microstructure. There is no physical length scale assigned to the voxel edge length, since these models are intended to be general. In real composite materials, the size of inclusions can range from tens of nanometers (gold nanoparticles in polymer matrices) to tens of millimeters (gravel in concrete). Digital spheres of diameter 41 voxels were used. They were defined by centering a continuum sphere of 41 voxel diameter on a voxel and labeling each voxel whose center fell in the interior of the continuum sphere as belonging to the sphere. The voxel volume turned out to be 36,137 voxels, which only differs by -0.1% from the continuum value of 36,086.95. A list of voxel centers was made, in reference to the center of the sphere, which was used to insert a sphere into the unit cell. Spheres were inserted in two ways. The first way was inserting the spheres, centered at randomly chosen voxel centers, but not allowing them to overlap. Periodic boundary conditions were used so that a part of a sphere that extended beyond the unit cell boundary was filled in on the other side of the unit cell. Figure 2a shows an 800 9 800 pixel slice of this microstructure, called the two-phase non-overlapping sphere (2NOS) model. Note that the ratio between the size of the unit cell and the size of the inclusions, about 20, means that finite size errors will be negligibly small, since a ratio of 5 or 10 would be adequate for having negligibly small finite size errors [21-24, 29, 30]. Using a ratio of 20 assures us that finite size errors can be completely ignored. The circles in Fig. 2a show a size distribution since they are sliced through at different distances from monosize sphere centers. This model has a two-phase structure, consisting of the percolated matrix (phase 1) and the isolated volume inside the spherical inclusions (phase 2). The volume fraction of phase 2 is then just the number of spheres, N = 4719, times the volume per sphere (36 137 voxels), divided by the unit cell volume. In the case studied, the sphere volume fraction was 0.333 and the matrix phase fraction was 0.667. The second method of constructing the model microstructure was to insert the spheres, centered at randomly chosen voxel centers, but allowing them to freely overlap. In a two-phase infinite system of matrix plus identical overlapping objects, randomly centered, the remaining matrix volume fraction is given by [31][32][33]: where V o = the volume of one object and n = the number of objects per unit total volume. For any finite realization of this kind of microstructure, even if continuum rather than digital objects are used, Eq. (3) will be an estimate only. The resulting two-phase microstructure, called the two-phase overlapping sphere model (2POS), had phase fractions similar to the 2NOS model. A slice is shown in Fig. 2b. 5738 spheres were placed, so that Eq. (3) gives a prediction of 0.6670, which was very close to the actual matrix phase fraction of 0.6656. A burning algorithm was used to determine that both phases were percolated [34]. The second phase was expected to be percolated since its volume fraction, 0.3344, was greater than 0.2895 ± 0.0005, the percolation threshold determined for the 2OS model [35].
To build a simple three-phase overlapping sphere composite model (3POS), the overlapping sphere algorithm was slightly modified. When each sphere was randomly inserted into the matrix (1), the sphere (meaning all volume it replaced) was chosen to be phase 2 with a probability of 0.3 and phase 3 with a probability of 0.7. The final volume fractions actually achieved were: / 1 = 0.3500, / 2 = 0.1962, and / 3 = 0.4538. In the system considered, nV o = 14 850 V o /(800 3 ) = 1.0481142, which gives a matrix phase fraction of 0.3506, which differs by less than 0.2% from the achieved value, probably due mainly to the small digital error in using digitized spheres vs. continuum spheres [36] and the fact that the model is only a finite-size realization. The total inclusion volume fraction is predicted by Eq. (3) to be 0.6494, which is very close to the sum of the volume fractions of phases 2 and 3 (0.6500). In the 3POS model, the phase 2 and 3 volume fractions are predicted to be 0.6494 Á 0.3 = 0.1948 and 0.6494 Á 0.7 = 0.4546, which are again very close to what was actually found. Since the matrix phase fraction, 0.35, in the 3POS model is much less than the matrix phase fraction, 0.667, in the other two-phase models, the inclusion phases, which have higher moduli than does the matrix, will carry more weight in the average E value. A burning algorithm [34] showed that phases 1 and 3 were percolated and phase 2 was not. A percolation phase diagram for nPOS models, where n C 3, is discussed in more detail in Ref. [37]. Figure 2c shows one of the possible 800 9 800 voxel cross-sections of the 3POS model. In the first columns of Table 1, the elastic moduli assigned to each phase and the phase volume fractions in all three models are given. Note that there is perfect contact between voxels of the same or different phases. The last column, marked ''Unit cell,'' shows the effective moduli obtained by the FEM computation averaged over the entire unit cell. Note that in experiments on multiphase materials, a given phase may itself have a distribution of elastic moduli, centered around an average value, due inherent randomness in an amorphous phase or crystallite orientation in a polycrystalline phase. In this paper, each phase has a distinct value of elastic moduli. It will be seen below that the E distribution histograms arise only from mixing of these precise phase moduli values, not from any distribution of elastic moduli values inside a given phase. It is interesting to note that the computed composite values of E and m, listed in Table 1, agree rather closely with the self-consistent effective medium theory (SC-EMT) listed on p. 475 of [1]. The SC-EMT values for both two-phase models are E = 16.77 GPa and m = 0.267, and the three-phase model has SC-EMT values of E = 39.63 GPa and m = 0.226.
To better understand how the size of the SVs affects the E distributions, one needs a length scale inherent to the phase in the three models in order to understand how the E distributions change with SV size. Comparing Fig. 2a, b, it is reasonable that the size scale of the inclusion phases in the 2POS model could be somewhat larger than in the 2NOS model, since overlapping spheres contain larger lengths in them than do isolated spheres. However, there is no obvious intuition for the 3POS model, as seen in Fig. 2a, c common length scale is thus needed to be able to quantitatively compare the three models to each other. Reference [1] gives a choice of several microstructuredependent length scales, which could be used for our isotropic microstructures, but none are particularly well-suited to the oriented cubic SVs used in this paper. The closest length-scale-type quantity is the coarseness, defining how the phase volume fractions fluctuate over local length scales [1,38,39], but is still not entirely appropriate since our SVs are aligned with the voxel lattice. In this paper, and in analogy to the lineal path function [1], we define the quantity W j (q) to be the probability that a randomly placed q 3 voxel cube, oriented with the 3D voxel arrangement, will be entirely contained in the jth phase. Since the limit q ? 1 is a single voxel, the value of W j (1) is just / j , the volume fraction of the jth phase, since in this limit W j (q) just adds up the number of voxels belonging to phase j. For some value of q, called q j0 ,W j (q) will go to zero and remain zero for larger values of q, since there would be no place in the microstructure where this size cube or larger could fit anywhere into the jth phase. The average value of q for phase j, computed using W j (q) and Eq. (4), will then be a length that is representative of phase j, where the largest value of q considered was Q. However, as will be seen below, the q j0 parameters are actually more useful than the averages for understanding how the E histograms change.
A simple way of computing W j (q) is to randomly overlap an oriented q 3 voxel cube on top of the microstructure and sum the values of the voxels overlaid (1 = phase 1, 2 = phase 2, 3 = phase 3). If this sum is equal to q 3 , then the cube is entirely in phase 1; if the sum is equal to 2q 3 , the cube is entirely in phase 2; if the sum is equal to 3q 3 , the cube is entirely in phase 3. One million random points were chosen for each model microstructure, which was Table 1 A list of the elastic moduli (units of GPa) assigned to the various phases (1, 2, or 3) in the three microstructures investigated, along with the phase volume fractions and the values of the average over the entire unit cell determined to be sufficient by comparing to W j (q) computed by systematically summing over every voxel in the 800 3 models. An important number for each model will be q = q j0 , which is that value of q such that W j (q) is less than 0.0001, a reasonable approximation of where W j (q) is zero considering the finite size of the unit cell of the models, which would tend to make the large q tail of W j (q) appear longer than it really is. If one were using some kind of microstructural probe, such as a nanoindenter, then q j0 gives the largest size at which the probe might randomly fall inside a single phase.

Results
Each microstructural system was solved for the nodal displacements, and the six stresses and six strains were averaged over each voxel. into the given E bin. Not all the graphs corresponding to all the choices of SV sizes will be presented, but only enough so that the changes in the E histograms can be followed as SV size increased.

W j (q) versus q
We first present the W j (q) results, which will be used to understand the E histograms. Figure 3 shows plots of W j (q) for (a) the 2NOS and 2POS models, and (b) the 3POS model. Notice that in Fig. 3a, the W j (q) curves for the 2POS model are always larger than those for the 2NOS models. Table 2 lists the values of \q[ and q j0 computed for the different phases and models using the W j (q) curves from Fig. 3.
In Table 2, and referencing Fig. 3, note that that phase 1 in the 2POS model has significantly more ''size'' than does phase 1 in the 2NOS model, while the sizes of phase 2 in both models is quite similar, in terms of \q[ and q 0 . The ''size'' of phase 1 in the 3POS model is smaller than both the 2NOS and 2POS models, and phases 2 and 3 are similar to the twophase models.

Two-phase non-overlapping sphere (2NOS) model
As given in Table 1, the global effective E was found to be 16.30 GPa. This number will affect how the following E histograms transform as the SV size increases, since the distribution goes to a delta function of weight 1 centered at 16.30 GPa in the limit where the SV becomes equal to the size of the unit cell. Figure 4a shows the E distribution computed over 4 3 voxel SVs. There are clear, sharp peaks at values of 10 GPa and 50 GPa, which are the E values of phases 1 and 2, respectively. The values of the phase fractions at these peaks show that they are made up of a single phase only. However, there is a small amount of weight for values between these two peaks, and the phase fraction curves show how the phases are starting to mix to give the intermediate E results. For example, the bin at 24 GPa is made up of almost exactly 50% of phase 1 and 50% of phase 2. If a sampling instrument had linear dimension of 4 voxels or less, and many probes were made to give statistical results, the data would readily identify two major phases. Figure 4b shows the E distribution at the SV size of 16 3 voxels. The heights of both sharp peaks have decreased by about the same factor, while the weights of the intermediate E values, especially close to the 10 GPa peak, have significantly increased. The phase fraction curves show the phase mixing that is giving rise to this additional weight near 10 GPa.
The next two graphs in Fig. 4 have SV sizes that are closer together, since at this point qualitative changes take place rapidly in the histograms as a function of SV size. Figure 4c, for a SV size of 24 3 voxels, shows that a small peak still remains at 10 GPa, but the peak at 50 GPa has mainly disappeared. There is still some weight there, since there is still a bin at 50 GPa and the phase 2 fraction still goes to 1 as E goes to 50 GPa. However, at least 80% of the weight of the histogram is in the wide shoulder between 10 and 35 GPa. The phase fraction curves look quite similar to those in Fig. 4a, b, although the actual phase volumes associated with each symbol are quite different, since they reflect the number of SVs with local E values that fall into a given bin of the E histogram. At this point, if a sensing instrument sampled a region with linear dimension equal to 24 voxels, it would not identify two major phases. Figure 4d shows that by a SV size of 28 3 voxels, there are no more pure phase 2 SVs, since the upper limit of the graph is at approximately 48 GPa and the phase fraction points for phase 2 now no longer go to 1 at the right hand side of the graph. This is predicted by the value of q 20 in Table 2, 24 voxels. There is no way to get a local SV E value of 50 GPa by any mixing of phase 1 and phase 2. The peak at 10 GPa has broadened and now values of up to 15 GPa have equal weight as the 10 GPa peak. There are still pure phase 1 bins, since the value of q 10 is 34 voxels and 28 \ 34.
Figures 4e, f show the E histogram at SV sizes of 32 3 voxels and 40 3 voxels. In Fig. 4e, the histogram has changed qualitatively and is now a single, asymmetric peak centered at around 14 GPa. There is still weight at 10 GPa, and the phase fraction for phase 1 shows that there is still a small amount of pure phase 1 SVs. The right-hand limit on the E axis has decreased to slightly less than 40 GPa, and there are less phase 2 and more phase 1 voxels at this limit.   Fig. 3 The computed W j (q) curves for (a) the 2NOS and 2POS models, and (b) the 3POS model   The global effective E, averaged over the entire unit cell, was listed in Table 1 as 16.69 GPa, about 2.5% larger than the 2NOS case, so the center point of the eventual single peak in the E distribution will be close to that of the 2NOS sphere model. Figure 5a shows the distribution of E computed over 4 3 voxel SVs. There are clear, sharp peaks at values of 10 GPa and 50 GPa, which are the E values of phases 1 and 2, respectively. The values of the phase fractions at these peaks show that the peaks are made up totally of a single phase. However, there is a small amount of weight for values between these two peaks, and the phase fraction curves show how the phases are starting to mix. Figure 5a shows that if a sampling instrument had linear dimension of 4 voxels or less, it would readily identify two major phases. Figure 5b skips to the 32 3 SV results. The SV sizes between 4 3 and 32 3 voxels show a progression that is quite similar to that shown for the 2NOS model, so they will not be presented. For the 32 3 SV results, the 50 GPa peak has vanished, since q 20 = 26 voxels for the 2POS model and 32 [ 26. There still could be a small amount of weight showing in the last histogram bin, however, because the bins are of finite size. The peak at 10 GPa is still visible, although much reduced in height, and there is a large shoulder extending to about 40 GPa. Figure 5c shows the E distribution for a slightly larger SV size, 36 3 voxels, still showing a distinct phase 1 peak. Figure 5d, for a 40 3 voxel SV, shows that the 10 GPa peak is no longer a distinct peak. There is now a single peak, gradually decreasing from 10 to 40 GPa, but there is still weight in the 10 GPa bin, since q 10 = 56 for this model and 36 \ 56.
At an SV size of 80 3 voxels, as shown in Fig. 5e, there is a clearly defined single peak, centered at about 17 GPa, and the phase fraction lines have become much flatter, in terms of their minimum and maximum values, than in Fig. 5d., and there is no weight at E = 10 GPa, as expected, since 80 [ q 10 = 56. A sensing instrument of linear size greater than about 40 voxels, which is about the size of one inclusion, would only see one peak and therefore one phase but with a spread of E values. A sensing instrument with linear size larger than that, up to 80 voxels in size, would show a single phase, with a distribution of elastic moduli values, which was centered on the global composite modulus value.

Three phase overlapping spheres (3POS) model
This model was chosen so that more than two phases, with a wider range of phase elastic moduli, could be studied. More theoretical work can be done with only two phases [1], but many composite materials have more than two phases. In Table 1, the global effective E, averaged over the entire unit cell, was found to be 37.56 GPa, which is substantially greater than the twophase models. Figure 6a, b show the E distributions and accompanying phase fractions for SV sizes of 4 3 voxels and 16 3 voxels. In Fig. 6a, three clear, sharp peaks are seen at 10 GPa (phase 1), 60 GPa (phase 3), and 150 GPa (phase 2), corresponding to the phase E values, and only one phase fraction is non-zero at each peak. The phase fraction curves for phases 1 and 2 both have cusps at the 60 GPa peak, so that only phase 3 is nonzero there. In Fig. 6b, for a 16 3 voxel SV, the peaks have shrunk in height by factors of 2 to 3, with the greatest relative shrinkage being in the middle at the 60 GPa peak. Significant weight has developed between the 10 GPa and 60 GPa peaks and between the 60 GPa and 150 GPa peaks, showing the phase mixing that has taken place with the increased SV size.
For a SV of 24 3 voxels, there has been a sharp qualitative change in the E histogram, as can be seen in Fig. 8c. The 150 GPa phase 2 peak and any weight at all in the last bin has essentially disappeared, since q 20 = 24, and a distinct 10 GPa phase 1 peak has almost disappeared, too, though with weight still in the 10 GPa bin. The phase 3 peak can still be seen at 60 GPa, but it is no longer pure phase 3. Most of the weight of the distribution is between 10 and 60 GPa, since the eventual limit of the distribution is a narrow peak at 37.56 GPa. In Fig. 6d, at a slightly larger SV,  3 voxels, there is now a single broad peak, and the phase 3 bin at 60 GPa is even more mixed phase. Note that the 60 GPa peak in Fig. 6c is only about 70% phase 3, but it still exists since there are many ways 150 GPa and 10 GPa voxels could mix to produce an E value of 60 GPa. In the 60 GPa peak, we have contributions to the ''real'' phase 3 peak by combinations of phases 1 and 2. This means that even if peaks can be recognized in a technique like statistical indentation, since they are appearing at E values that are close to those expected for a pure phase (so they can be labeled), no volume fraction of the phases can be calculated from these peaks. Taking a closer look at the data shown in Fig. 6c, the value of the phase 3 volume fraction was averaged for all the 24 3 voxel SVs that had a value of E between 50 and 70 GPa, in 2 GPa-wide bins in that range, with the results shown in Fig. 7 plotted against the mid-point E value of each bin. The open circles show the maximum value of the phase 3 volume fraction encountered over the bins examined, some of which were near unity, so that some of the SVs in this region were indeed almost pure phase 3, which is reasonable, since 24 \ q 30 = 27 voxels. However, for every bin in the histogram, the minimum value of the phase 3 volume fraction in an SV was 0 (not shown in figure), indicating that there were some SVs with values of E in this range that contained no phase 3 voxels. The open squares show the average value of the phase 3 volume fraction in each E range, which peaks near E = 60 GPa as might be expected. The solid circles show the number fraction of the SVs that had a phase 3 volume fraction  Figure 7 shows that there was a significant number of SVs included in the peak around 60 GPa whose local microstructure had very little to do with phase 3. Counting these contributions to the peak at 60 GPa as indicating the presence of phase 3 would therefore be erroneous. Note that the 60 GPa peak actually became less than 100% phase 3 by about an SV size of 20 3 (data not shown). By the time a 32 3 voxel SV is reached, Fig. 8a shows that the single peak has become better defined but has about the same width as in the 28 3

voxel SV.
There is still some weight in the 10 GPa bin, but this could be just because the bin width in Fig. 8a is larger than in Fig. 6, since q 10 = 31 voxels, which is close to 32. In Fig. 8b, using a 40 3 voxel SV makes the single peak significantly narrower and somewhat higher. The phase fraction curves have become much flatter, too, and show a large amount of mixing over the entire E range. There is very clearly now no pure phase 1 bin, since 40 [ q 10 = 31. A sensing instrument of linear size greater than about 32 voxels, which is about the size of one inclusion, would only see one peak and therefore one phase but with a spread of E values not centered on the global composite value. A sensing instrument with linear size larger than that would show a single phase with a steadily narrower distribution of elastic moduli values, which was more and more centered on the global composite modulus value.

Discussion and future work
Many figures have been shown in this paper of how the morphology of the E distributions change with the SV length scale. It should again be noted that these E distributions come about only because of microstructural mixing within a given size SV, since the phase moduli assigned to each phase was exact and was not distributed within a phase [9].
Although the elastic fields around any kind of probe will differ, in general, from simply summing the local stress and strain fields in a SV of a material to which a uniform strain was applied, the results obtained are still instructive in that the development of the E histograms with SV size should probably be qualitatively similar. For these kinds of models, the function W j (q), which is the probability of finding an oriented q 3 cube entirely in phase I, and the values of q j0 , the value of q at which W j (q) approximately goes to zero for phase j, can be used to understand how the E histograms change with SV size. The value of q j0 predicts where peaks associated with phase j drop out of the E histogram for the two-phase models, since there was no combination of phases that could produce the E values associated with the phases. However, it was shown for the 3POS model that a distinct peak, located at the phase 3 (E = 60 GPa) bin, could reflect a mixture of phases and not a pure phase 3 peak, as determined by numerical simulation and the value of q 30 . This is because a mixture of 10 GPa and 150 GPa material could produce a 60 GPa effective E in a given SV. For each microstructure, the approximate size of the probe was found such that only a single peak would be seen, with its center approximating the effective E of the entire system. The width of the single peak does not reflect anything about the actual number of distinct phases in the system [9].
It should be noted that the models and data analysis presented in this paper were ideal, in several ways, compared to an experimental approach. First, the nanoindentation and atomic force microscope shown. The number fraction of the SVs that had a phase 3 volume fraction less than 20% in each E range (solid circles), and the number fraction of the SVs that had a phase 3 volume fraction greater than 80% (solid squares) are also shown approaches only sample the surface of a material, so that much less of the material is sampled than in these models, where the entire volume was used. A consequence of this is that much fewer SVs are used in experiment to determine E histograms than in these models. Of course, larger areas could be examined to partially make up for this surface/volume difference. Note that if there are inclusions in the material, then in general, a cut through the material to produce a smooth surface will produce an apparent size distribution of inclusions that range down to smaller size inclusions [40,41]. This would probably tend to make the SV size decrease that was needed to make the multi-peak structure disappear [42]. Second, determining the local moduli was done exactly, albeit numerically, in these models, while in experiment the local moduli have to be extracted using assumptions such as an assumed value of local Poisson's ratio [6] on a surface that may have been damaged by polishing and sample preparation. Indeed, a substantial fraction of measurements typically needs to be discarded because of damage on the surface (e.g. [17]). Roughly speaking, the results of these models can be stated that by an SV size of about one half the diameter (41 voxels) of the spherical building block of the models, one could not identify more than one peak from any of the models.
However, such small indents may be impossible on some composite materials such as cement paste, since the intrinsic surface roughness of the paste might be of the same size and cannot be polished out in a nanoindentation experiment [43]. Third, the local ''probe'' used in this paper was ideal, while that used in experiment has uncertainties associated with it in terms of quantities such as assumed shape and calculated penetration depth [6][7][8]. Fourth, while in the model the phase make-up of every SV is known, in mechanical measurements on the surface the local phase make-up, especially below the visible surface, is almost never known. In experiment, even if a phase can be clearly identified on the surface, the probe may be sampling other phases just below the surface of the sample [44]. In recent practical applications of statistical nanoindentation, researchers have started to avoid an approach based only on the analysis of the mechanical response. When extra microscopic information is available about the phases present within the interaction volume of each indentation point, it is possible to discard points where multiple phases are present and finally identify the mechanical response of the single phases [45]. Of course, discarding the majority of points (either because they correspond to a phase mixture or because the surface is damaged) requires a larger number of indentations and precludes the possibility of deriving the phase fraction of the different phases.
The type of materials that the analysis of this paper could apply include any multi-phase composite material. Examples would include cement paste [8,42,43] or concrete [46], or a polycrystalline, poly-phase  [47,48] or other alloys. Polymers filled with an inorganic powder would also lend themselves to this kind of analysis. Future work will look at real cement paste microstructure in 3D, as captured by focused ion beam serial sectioning [42,43], and how variations in the microstructure, elastic moduli contrast, and data analysis affects the clarity of the Young's modulus and phase volume fraction histograms. Specifically, three areas will be studied: the effect of changing from monosize inclusions to an inclusion size distribution, the effect of varying the elastic moduli contrast between phases, and the effect of having ''coating'' phases of a known thickness [49].
Funding Open Access funding provided by Lib4RI -Library for the Research Institutes within the ETH Domain: Eawag, Empa, PSI & WSL.

Compliance with ethical standards
Conflict of interest The authors declare that there is no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.