Effective hyperelastic material parameters from microstructures constructed using the planar Boolean model

We construct two-dimensional, two-phase random heterogeneous microstructures by stochastic simulation using the planar Boolean model, which is a random collection of overlapping grains. The structures obtained are discretized using finite elements. A heterogeneous Neo-Hooke law is assumed for the phases of the microstructure, and tension tests are simulated for ensembles of microstructure samples. We determine effective material parameters, i.e., the effective Lamé moduli λ∗\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda ^*$$\end{document} and μ∗\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu ^*$$\end{document}, on the macroscale by fitting a macroscopic material model to the microscopic stress data, using stress averaging over many microstructure samples. The effective parameters λ∗\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda ^*$$\end{document} and μ∗\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu ^*$$\end{document} are considered as functions of the microscale material parameters and the geometric parameters of the Boolean model including the grain shape. We also consider the size of the Representative Volume Element (RVE) given a precision and an ensemble size. We use structured and unstructured meshes and also provide a comparison with the FE2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^2$$\end{document} method.


Introduction
The problem of relating effective constitutive parameters of composite material laws to their microstructure is common in many areas, e.g., in material science as well as in geoscience. Approaches to calculate the elastic moduli of aggregates reach back over a century; see [20] and the references therein. In the regime of linear elastic behavior many well-known contributions have been made to predict effective elastic moduli of multiphase materials using bounds [2,17,18,52] or estimates [22,65]; also cf. [16].
It is well known that effective properties of composite materials do not only depend on the mechanical properties of the phases and their volume fractions, but also on other details of the microstructure, e.g., the shape of the inclusions, their orientation, the connectivity of the phases; cf. [16,34,80,90].
In this paper, 1 we consider two-phase random microstructures, constructed by stochastic simulation using the Boolean model [8], which is a random set model, which models inclusions in a matrix material. These grains may have different shapes and sizes and also their density is a model parameter. Throughout this paper, we always consider stiff inclusions embedded in a softer matrix. The microstructure samples are discretized using finite elements and we use a nonlinear Neo-Hooke material model. By applying numerical homogenization, using many microstructure samples, we obtain effective material parameters on the macroscale for various parameter sets of the Boolean model. This enables a systematic study of the effective material parameters, i.e., of the Lamé constants λ * and μ * , as a function of the parameters of the Boolean model.
We consider different grain types, i.e., discs and rectangular inclusions with different aspect ratios. The numerical homogenization approach is based on applying parameter identification to the microscopic problem, assuming a given material model on the macroscale. The parameters on the macroscale are identified by minimizing the squares of the residuals of the equilibrium equation, which is also known as equilibrium gap method (EGM) [9]. A sufficient fit could be achieved assuming a Neo-Hooke material model on the macroscale. This has been observed earlier by Löhnert et al., who have also fitted a macroscopic Neo-Hooke law to microscopic data [58][59][60]. Note that similar fitting approaches have also been suggested for much more complicated material models such as crystal plasticity [54,55].
We will also use identified macro material parameters to perform a comparison with simulations using the FE 2 computational homogenization method. In the FE 2 no macroscopic material law needs to be defined. The tensors describing the effective properties of the microstructure are directly coupled to the macroscale; see, e.g., [76]. The FE 2 method is computationally expensive but parallel scalability to more than a million parallel processes has been demonstrated [47]. Also FE-FFT based two-scale methods have been developed, where the microscopic problem is solved using fast Fourier transform; see, e.g., [48,75]. For a comprehensive view on homogenization in solid mechanics, we refer to [43,64,73,76] and the references therein.

Related work
Our paper follows the example of Jeulin [37] and his coauthors in using models of stochastic geometry such as the Boolean model for the construction of random microstructures by stochastic simulation and investigation of their material properties, [10,[30][31][32][33][34]40,42]. The Boolean model has been used even earlier in the context of porous media by Matheron; see [77] and the references therein. It was subject of studies in the context of linear elasticity [29,85], permeability [84], and conductivity [72]. Jeulin and Willot have worked on conductivity and nonlinear elasticity using a power-law; cf. [86]. In [46], Neo-Hooke and Ogden-type material models with a contrast, i.e., a ratio between the parameters of the phases of, up to 10 were considered in the context of estimating the size of the representative volume element. The authors studied a two-dimensional random composite of randomly distributed non-overlapping discs of constant diameters. They also used a fitting procedure to determine macroscopic material parameters for Neo-Hooke or Ogden material laws. In [58,59], non-overlapping spherical random inclusions were considered using a Neo-Hooke law, also using a contrast in the material parameters of 10.
Many other authors have considered overlapping inclusions and we also follow this approach by using the Boolean model. This will lead, for a fixed area fraction, to a stiffer macroscopic material response since large stiffening structures can emerge. Of course, overlapping grains are not always realistic, as for example in concrete.
The prediction of the effective material behavior of composite materials has been studied extensively during the last century. Analytical homogenization approaches for linear composites go back to the works of Eshelby [13], Kröner [51] and Hill [22]. Corresponding approaches to the determination of limits and approximations for anisotropic properties and constituents were presented by Willis [82,83] and Hashin [15]. In the theory of nonlinear composites Castañeda et al. [3][4][5]7,56,57] and Suquet et al. [5,62,78] have contributed variational procedures, improved bounds and homogenization methods for the prediction of effective nonlinear material behavior. For the special case of strongly anisotropic phases Idiart et al. [28] and Willot et al. [87,88] extended these methods and investigations.
In homogenization the notion of a representative volume element (RVE) [21,37] is central. The RVE is required to statistically represent the microstructure. Determining the correct size of the RVE has been a well-studied topic of the effective continuum theory; see, e.g., [26,39,42,68,79] in the linear elastic and [14,45,46,70,79] in the nonlinear elastic regime. We use the approach described in [42] to determine the RVE.
For problems where the scales become comparable, second-order homogenization methods have been introduced; see [50,73] for a comprehensive view on the topic. Analytical homogenization methods are mainly used in the linear elastic case; cf. [43] and the references therein. For nonlinear composites Castañeda and Suquet developed ana- Fig. 1 In order to illustrate the discretization, we show both the random structure and the finite element mesh. The figure displays two small microstructure samples (ratio of volume element size and grain size is small) of Boolean models with circular (top) and rectangular grains (bottom). Both samples belong to an area fraction of A A = 0.5. The left images display grains (blue) and germs (red). The resulting pixel data is shown in the middle. The corresponding finite element mesh, using four linear elements for each pixel, is shown on the right. The volume elements used in our computations are significantly larger, e.g., by a factor of 256; see lytical homogenization methods in the hyperelastic case; see [43,73] for an overview on that topic. In [6], Castañeda and Tiberio introduced the idea of a linear elastic, homogeneous comparison composite to estimate the effective constitutive behavior of the hyperelastic composite with heterogeneous microstructure; see also [24].

The Boolean model
In this paper, we use the Boolean model [8, pp. 64 ff] to construct microstructures. The Boolean model is a standard model of stochastic geometry introduced by Matheron; see, e.g., [8,63] and the references therein. Geometries constructed by the Boolean model can model very disperse spatial patterns and are suitable for the simulation of interconnected composite materials due to the overlapping of its grains; see [30][31][32][33][34]40,77]. The Boolean model has many application in physics, e.g., in permeability problems [84], in conductivity and porous media [72], and in material science [29,35,36,41,61,85,86]; see also [8] for an overview. A Boolean model is a random set in R d and as such suitable for modeling bi-phase structures, where the set constitutes phase 1 and the complement phase 2. We work here in the planar, two-dimensional case.
The Boolean model Ξ is a special germ-grain model and is constructed by germs x n , which form a homogeneous Poisson point process, and grains, which are compact sets Ξ n . The Ξ n are independent identically distributed compact subsets of R 2 ; a further set with the same distribution is called the typical grain and is denoted as Ξ 0 . In this paper, the grains form the stiff inclusions embedded in a softer matrix.
The corresponding Boolean model is the random set Note that the sets Ξ i + x i may overlap, see Fig. 1. The set Ξ is statistically homogeneous, i.e., its distribution is translation invariant. In the following, we assume that Ξ 0 has a rotation-invariant distribution. Then Ξ is also isotropic.
We denote the intensity of the Poisson process by θ , which is the mean number of germs per unit area. As grain types we use discs, denoted by D, and rectangles, denoted by R 1 and R 3 , with aspect ratios 1 : r r , where r r = {1π, 3π }. Microstructure samples of the planar Boolean model using different grain types are shown in Fig. 2.
Two important global characteristics of a Boolean model are the area fraction A A of the inclusions (the area fraction of the matrix phase is (1 − A A )) and the specific boundary length L A [8, p. 80]. The area fraction A A is the fraction of the plane covered by Ξ or the probability that the origin belongs to Ξ . It holds where E is the expected value and ν 2 is the planar Lebesgue measure. The specific boundary length L A is the mean length of the boundary ∂Ξ of Ξ per unit area. It holds Here, E(ν 1 (Ξ 0 )) is the mean boundary length of the typical grain and ν 1 is the Lebesgue measure.

Constructing microstructure samples using the Boolean model
In this paper, the microstructure samples of the Boolean model are always constructed on = [0, 1] 2 . Instead of constructing a larger microstructure sample on a larger domain we consider smaller grains and a higher intensity θ of the Boolean model. We determine the model parameters for the construction of the sample by, first, choosing a size of the microstructure sample and computing the size of the rescaled grains on = [0, 1] 2 , and second, given an area fraction of the Boolean model, we determine the intensity θ using (2). Our microstructure samples are periodic, i.e., grains overlapping the boundary of are continued on the opposed side. In [81], the effect of grains crossing or not crossing the boundary is investigated. The results show the importance of allowing grains to cross the boundary of .
For each set of microscale model parameters, we construct n microstructure samples, which we refer to as microstructure ensemble. We use the terms microstructure sample and volume element (VE) as synonyms.
Except where mentioned explicitly, we process the microstructure samples of the Boolean model as black (matrix) and white (inclusion) pixel data. The pixel data is then discretized using linear finite elements; see Fig. 1. Note that the low-order finite element discretization introduces significant discretization errors. However, higher order elements may not help in our regime with high contrast in the parameters. Instead, adaptive mesh refinement could be applied.
To construct the pixel data, we start with an empty image of N 2 pixels. Then, if the center of a pixel is overlapped by a grain, the value of this pixel is set to one. Figure 1 displays examples of microstructure samples of planar Boolean models and their representations as pixel data. Where we use unstructured meshes, they are constructed using the Matlab mesh generator Mesh2D [11].

Homogenization using parameter identification
In our computational homogenization approach, we use a hyperelastic, homogeneous comparison material and identify its parameters by minimizing the squares of the residuals in the equilibrium equation. This parameter identification method is well known in mechanics and often denoted equilibrium gap method (EGM) [9]; see Sect. 4.2. However, in our homogenization procedure, we use the mean stress, averaged over the microstructure ensemble, as reference in the minimization; see also [58][59][60]. An advantage of this approach is that the (expensive) parameter identification is carried out not for every microstructure sample but only once for each ensemble of n microstructure samples.

Load cases for the parameter identification
For our parameter identification, we consider uniaxial tension tests with 10% deformation applied to the microstructure samples; we then consider a biaxial test and a shear test for the verification of the parameter fit; see Fig. 3. It is well known that the choice of the boundary conditions has a considerable influence on the effective material tensor fields of a composite and, if the volume element is too small, introduces a significant bias; see also Sect. 4.4. In the elastic context, this is related to the Hill condition or the Hill-Mandel condition see, e.g., [19,21]. To reduce the effect Ω Ω Ω Fig. 3 The boundary conditions for the uniaxial tension in x-direction (left), the non-proportional biaxial tension (center) and the shear test (right). Dirichlet boundary conditions are displayed by triangles. On the remainder of the boundary, we apply homogeneous Neumann boundary conditions of the (artificial boundary) conditions, we need to consider volume elements of sufficient size (representative volume element). For a detailed computational investigation of this matter in finite elasticity, we refer, e.g., to the work of Ostoja-Starzewski; cf. [46,70]. In [46] the authors investigate the effective behavior of random two-phased composites using (different) hyperelastic material models for the phases. They found computationally for their hyperelastic models that Dirichlet and Neumann boundary conditions provide upper and lower bounds to the effective stress and the strain energy. They also show that orthogonal-mixed boundary conditions yield an intermediate result for the effective material behavior.

Parameter identification using least squares for the averaged stress fields
Throughout this paper, we use an isotropic Neo-Hooke strain energy function as the material law on the microscale, see, e.g., [92] and [91], The parameters λ and μ are the Lamé constants [25] on the microscale. We refer to the material parameters on the macroscale obtained by parameter identification as the apparent parameters, and we denoted them by λ and μ.
If the RVE criterion is satisfied (see Sect. 4.4) we refer to them as effective material parameters rather than apparent parameters and denote them by λ * and μ * .
Here, J = det(F), where F is the deformation gradient, I 1 = tr(C) = tr(F T F) is the first invariant of the right Cauchy-Green tensor C. We consider Neo-Hooke because it is nonlinear but still simple. Other material laws could also be considered.
As material parameters for the matrix phase, we use λ = 100 GPa and μ = 75 GPa. For the inclusions, we consider material parameters which are larger by a factor c (for contrast), where c can be up to c = 100. Furthermore, we assume plane strain conditions. Figure 4 shows stress-strain diagrams for the uniaxial tension tests in x-direction, in y-direction, for a nonproportional biaxial tension test (15% in x-direction and 12% in y-direction) and a shear test (7.5%). Only in this figure, we consider 15% deformation instead of 10% deformation.
The averaged tensor field components P 11 , P 22 , and P 21 are shown as well as the ensemble average of the components P 11 , P 22 , and P 21 . The shear test and the biaxial test are provided as verification. The results in Fig. 4 show that, in our context, using Neo-Hooke on the macroscale leads to a sufficient fit, although small deviations are visible, which are larger for higher contrast c. Similar observations were already made in [46, p. 165] for c = 10 using nonoverlapping random inclusions and in the works of Löhnert et al. [58,59], where non-overlapping spherical random inclusions were considered. The parameter fit is performed using the ensemble average of the averaged stress fields, i.e., averaged over many (n ≥ 200) microstructure samples. Such methods, based on least squares fits of stress fields, are related to stress recovery approaches; cf., e.g., the overview paper [43] and the reference therein.
The method used here differs from the use mentioned in [43] in that our stress fields are averaged over the complete microstructure ensemble before the parameter fit. Note that the finite element simulations to compute the forward problem for the different microstructure samples are trivially parallel.
We compute the average stress P by integration over the microstructure samples, i.e., Averaged over the microstructure ensemble, the computed sets of macroscopic stress-strain data D k , D k = F k , P k , k = 1, . . . , M, with M as the number of load steps, are used to minimize the squares of the residuals [9] g(α) = k R(α) T R(α). (6) Here, R(α) is the vector of the differences between internal and external forces and α = (λ, μ) are the identified material parameters. It holds, Here, v i are test functions and T N is the traction vector. Solving the sequence of linear systems from the Newton scheme, the identified material parameters are obtained for each microstructure ensemble of the Boolean model. We use ensemble sizes n ≥ 200; see Table 2. Note that the α are unique if the (nonlinear) energy (α) is linear in the material parameters α; cf., e.g., [71], as is the case for Neo-Hooke. The parameters obtained from the identification are typically denoted as apparent (material) parameters. As a verification in Fig. 4 we compare the identified macroscopic material model also with the mean tensor field data from an independent nonproportional biaxial test and an independent shear test; see Fig. 3. The fit of the verification tests in Fig. 4 is not perfect, but, in our view, sufficient. To improve the quality of the fit, a more sophisticated material model on the macroscale could be used, e.g., by adding terms to the Neo-Hooke energy.
However, based on the results presented in Fig. 4, we decided to use the Neo-Hooke material model on the macroscale for all computations presented in this paper.

Choosing the size n of the microstructure ensemble
In our homogenization approach, we perform the parameter identification using average stress data computed over many microstructure samples, as in [58][59][60]. This is opposed to other approaches, where the parameter identification is performed for each sample, and where the parameters are then averaged, e.g., [42,53,85]. Note that the difference between averaging first, followed by the parameter fitting and, on the other side, fitting first, followed by averaging of the parameters is, in our experiments, very small. Indeed, for 200 microstructure samples of the Boolean model with discs and c = 100 and for an RVE the difference was smaller than the confidence interval by an order of magnitude. While the authors in [42,53,85] have considered the confidence interval for the parameters, we , instead, consider the confidence interval for the ensemble averaged stress data which is used to fit the macroscopic model: For the uniaxial tension test in x direction, we consider the confidence interval for the stress P 11 of the microstructure ensemble. Here, increasing the number of microstructure samples n reduces the width of the confidence interval.
We are interested in the confidence interval for the averaged stress data since this gives us an additional information to evaluate the parameter fits: If the confidence interval is Assuming a normal distribution of P 11 , the absolute and relative margin of error for our stress data P 11 from the uniaxial tension test are given by where z is the z-score of the normal distribution, μ P = P 11 is the mean and σ P = σ ( P 11 ) is the standard deviation. In this paper, we use z = 1.96 and the confidence level 95%. If we would assume, e.g., a log-normal distribution, then this confidence interval would be an approximation only, which for large ensembles, such as ours (n ≥ 200), is of sufficient precision [67]. In Fig. 5, the parameter fit of the macroscopic material to the stress data is shown, including the ensemble-averaged stresses and the confidence interval. We clearly see that the confidence interval is small compared to the mismatch of the fit. We can conclude that the mismatch of the fit can be expected to be the main source of error in our computations.
In our numerical computations, see Sects. 5.6 to 5.9, we choose a number n of samples for our microstructure ensemble heuristically, and we see, after having performed the computations, that the radii of the confidence intervals are within a desired range, i.e., have a radius smaller than 2.5%; see Table 2. This holds for rectangles well. For discs with log-normal radii, especially with a high variability in the distribution, this does not hold throughout our computations, see 5.10.2.

Choosing the representative volume element
In homogenization, typically, a representative volume element (RVE) is chosen. The notion of a RVE was introduced by Hill [21]: the RVE should be as small as possible and large enough to be representative of the material; see [14] for an overview on the topic.
The problem of choosing an appropriate size for the RVE has been studied extensively in the literature. For hyperelasticity, we mention the work of Ostoja-Starzewski et al. [70].
It is well known that the use of pure Dirichlet boundary conditions (also denoted kinematic uniform or KUBC) typically overestimates the stiffness of the microstructure. On the other hand, the use of pure Neumann conditions (also denoted static uniform or SUBC) typically underestimates the stiffness. Better approximations are usually obtained using periodic boundary conditions (also denoted PBC). In the FE 2 simulations in Sect. 5.1, we use periodic boundary conditions. In the uniaxial tension tests, which we use for the identification of the material parameters, we use Dirichlet on part of the boundary and homogeneous Neumann boundary conditions on other parts; see Fig. 3. This is valid if the microstructure samples are chosen large enough. For the comparison of FE 2 , making use of the same boundary conditions (i.e., periodic) would seem to be more canonical. We wanted, however, to use a consistent method in the EGM homogenization throughout the paper. If the microstructure sample is large enough, the influence of the boundary condition should be negligible.
For volume elements constructed from simulation of random media, we also have a random error, induced by fluctuations of the geometry of the microstructure. Both the systematic and the random error decrease for increasing volumes. For random media there is no single size of the RVE [69]. Instead, we expect to obtain material parameters with a certain precision if a given number of microstructure samples is combined with a certain size of the volume elements [42]. We choose a sequence of volume sizes and study the behavior of the apparent material parameters [8, p. 237].
It is convenient to use the non-dimensional size characteristic δ of a volume element (see [46]), defined by where d < L L macro , and where d is the size of the heterogeneity at the microscale, i.e., the diameter of the ball circumscribing the typical grain and L is the side length of the computational domain, which is the unit square in this paper.
We determine the RVE sizes using the approach presented in Kanit et al. [42]; see also Jeulin [37]. We use the microstructure ensemble data to identify the macroscale material parameters λ and μ for each microstructure sample and for different VE sizes to compute the variances D 2 λ (δ), D 2 μ (δ). Following [42], we then fit, first for the material parameter λ, the integral range δ 2 and the coefficient α to the data, using the model where d = 2. Then, we perform the same procedure for the apparent material parameter μ. The fit is presented in Fig. 18 in Sect. 5.5. A value of α close to 1 indicates a low bias in the microstructure ensemble data; see also [38].
The size of the RVE δ RVE can be defined for the parameter λ given a microstructure ensemble size n and a relative error Here, E(λ) is the mean of the apparent parameter ensemble. For the material parameter μ we define δ RVE (μ) and rel (μ) correspondingly. Throughout this paper, to define the size of the RVE for λ and μ, i.e., δ RVE (λ) and δ RVE (μ) for the approach in Kanit et al. [42], we will use rel = 0.01 and an ensemble size n = 200 as a reference. Indeed, n = 200 since this is the lowest number of realizations which we use; see Table 2. For clearness, we have refrained from presenting in our figures also other choices of rel and n for the approach in Kanit et al. [42] in addition.
Note that while we use Kanit et al. [42] to determine the RVE sizes, we compute our effective parameters as described in Sect. 4.2, i.e., we average the stresses before identification.
The integral ranges and RVE sizes are given in Tables 4  and 7 for different contrast and different grain types and in Tables 9 and 12 for discs with log-normal radii. In Table 8 we also show the RVE sizes for rel = 0.01 for different contrast and different grain types for n = 1.
The RVE sizes are also shown in our figures presenting the macroscale material parameters as functions of δ, i.e., Figs. 18, 21, and 25.
Plate with circular hole (left) and used symmetric part of the body and boundary conditions (right) If the RVE criterion is satisfied we speak of effective (material) parameters rather than apparent parameters. Here, we mean a low bias from the boundary conditions and a δ which is should be greater or equal to the maximum of δ RVE (λ) and δ RVE (μ) for rel = 0.01 and n = 200.

Numerical results
We use linear finite elements and use structured meshes (see Fig. 1), except where explicitly mentioned otherwise. We assume a Neo-Hooke material for the heterogeneous microstructures, using a contrast c, and identify the parameters of a Neo-Hooke material on the macroscale. Throughout this paper, the coefficient contrast c is identical for both Neo-Hooke material parameters λ and μ, i.e., the inclusions are stiffer by a factor of c. The random number generator of Mat-lab2019a is used. Note that for a high contrast of c = 100, we sometimes have a small number of microstructure samples where the finite element simulation failed to converge even for very small initial load steps in our adaptive load stepping scheme.

Comparison with the FE 2 computational homogenization method
Now, we compare the results of the homogenization procedure described in Sect. 4 with the FE 2 direct computational homogenization method, using a plate with a hole as a test case; see Fig. 6. The boundary value problem is first solved as a single scale computation using the material parameters obtained from homogenization. We then compare the results with the two-scale FE 2 scheme.
For this purpose a specific microstructure sample of the Boolean model is used; see Fig. 7. The material parameters of the two phases are chosen as λ = 100 and μ = 75 with a contrast of c = 100 for both microscale parameters. These parameters are then applied to a single scale computation of a quarter of the plate, discretized by 1543 linear triangular elements and the same Neo-Hooke energy function.
For comparison, the same boundary value problem is then solved with the FE 2 scheme as macro problem using the microstructure sample in Fig. 7 as underlying micro problem, constrained by periodic boundary conditions. The FE 2 computations have been performed using the AceGen and AceFEM packages (version 7.006), see [49], of Mathematica (version 12.0), see [89].
The (macroscopic) results of both computations are compared in terms of the von-Mises stress and the pressure as The corresponding contour plots are shown in Figs. 8 and 9.
The plots of both methods are in good agreement for the two quantities in characteristics as well as in amplitude even though the range of the chosen legend is quite small. This impression is confirmed by Fig. 10, is summarized in Table 1.
Despite this good correlation, it is notable that the FE 2 approach is advantageous in terms of localization analyses. The contour plots in Fig. 11 show the response of the micro structure corresponding to the specific integration point.  They reveal that the maximum von-Mises stress and pressure are significantly higher on the resolved RVE compared to the homogenized value on the macroscale. The gap is due to the combination of a high contrast in material parameters and a small volume fraction of stiff inclusions. In summary, if localization is important for the investigated application, it is reasonable to accept the considerably higher numerical costs of the FE 2 method.

Sanity check for the parameter identification
In Fig. 12, we consider, as a sanity check, volume elements using a centered disc, an off-centered disc, and a disc centered in the corner. We also consider the generalization using 64 discs. While we see that the results slightly differ for a single disc, for 64 discs the fitted stress-strain curves visually overlap.

Choosing the number of pixels per grain
For the first sets of computations, i.e., in the computations in Sects. 5.1 to 5.9, the microstructure samples of the Boolean model are given as pixel data; see Sect. 3.2. The use of such  In this paper, we use r N = 50 for most of our simulations data is standard in homogenization, especially if FFT-based solvers are used, e.g., [75,85]. We therefore first study the effect of the image resolution, i.e., we consider the apparent material parameters as a function of the grain resolution, r N , where E(ν 2 (Ξ 0 )) is the mean area of the typical grain and N 2 is the total number of image pixels. In Fig. 13, we show microstructure samples of Boolean models with discs and rectangles using different grain resolutions. Using the homogenization procedure described in Sect. 4, the apparent material parameters λ and μ are displayed for different values of r N in Fig. 14. For each data point n = 200 microstructure samples were constructed, and the apparent material parameters λ and μ obtained by the homogenization procedure are displayed as a functions of the grain resolution r N . We consider different values of the area fraction A A , two different contrasts, c = 10 and c = 100, and two types of grain, D and R 3 .
Note that the grain resolution directly translates to the size of our finite elements; see Fig. 1. In Fig. 14, we therefore observe the combined effect of a better resolution of the inclusions and a lower finite element error.
We also provide, as a reference, a comparison using unstructured meshes constructed using Mesh2D [11], as used in Fig. 15 (left), using dashed lines. Note that for Mesh2D, for discs of constant size, we use 30 points to approximate the boundary of the discs. This results in approximately 70 finite elements for each disk. For rectangles, we use at least 4 points for the shorter side.
From our results, as a compromise between accuracy and computational cost, we choose r N = 50 (see vertical line in Fig. 14) for all following investigations. For some cases with c = 100, e.g., discs with A A ≥ 0.5, a higher resolution would be desirable.
Using r N = 50, we tend to overestimate μ of our hyperelastic material model. For λ this does not hold in general. We find that for small A A and high ratio of L A /A A , we might tend to underestimate the apparent material parameters wrt. the grain resolution.

On the size of the representative volume element and on the size of the microstructure ensemble
In Fig. 16, the standard deviation of the apparent mean stress P 11 is displayed as a function of δ for different contrast as well as for different types of grain using high contrast.
Our results illustrate that ≈ Cδ −1 holds in our context for some constant C, i.e., for a hyperelastic material model using stiff inclusions, using different contrast and different types of grain; see Fig. 16. From (9) it subsequently follows that n ≈ Cδ −2 .

On the integral range and the size of the representative volume element
We consider the RVE size following the approach of Kanit et al. [42] (see also [37]) for discs with different contrast. For each microstructure ensemble we identify the ensemble of macroscale parameters and compute the variances D 2 λ (δ) and D 2 μ (δ) of both apparent material parameters λ and μ. We then use a least-squares fit (12) to the variances, identifying as parameters the integral range δ 2 and the exponent α; see  Note that, for the fit, we only use data points with small bias from the boundary conditions. Also note that, e.g., for discs we use ensemble sizes n as given by n = 200·(64/δ) √ 2 ; see Table 2 and its caption. Using sufficiently large ensemble sizes will result in small relative margins of error for the apparent material parameters and allow a good fit in (12).
In Fig. 18 we show the RVE sizes resulting from this approach for both apparent material parameters λ and μ and using rel = 0.01 and an ensemble size n = 200. In many of our computations, especially for small δ, we use ensemble sizes n which are significantly larger than n = 200, since our choice of the ensemble sizes n in Table 2 is independent of the area fraction A A . Using the approach of Kanit et al. [42] we see that a smaller number would be sufficient for certain values A A .
We see that for low contrast c = 2 and c = 10 using δ = 4 and δ = 16 is sufficient, even when using only n = 200. For high contrast c = 100, for two cases, i.e., for A A = 0.6 and A A = 0.7 for the material parameter λ a δ larger than δ = 64 would be needed according to this criterion. This is consistent with the confidence intervals for the stresses, which are significantly wider for c = 100.
Using this selection, we find α to be close to α ≈ 1 for all computations performed in this work except for discs with an extremely high variation in their radii; see Tables 4 and 9.
From these results, we can conclude that in sections 5.7, 5.9 and 5.10.2 our volume elements and ensemble sizes are chosen large enough and that we indeed compute effective parameters.

On the size of the representative volume element for discs for different contrasts
We would like to consider the size δ of the RVE, in the sense of Sect. 4.4, as a function of the parameters of the Boolean model. We use a grain resolution of r N = 50; see Sect. 5.3. For a given δ, the size of the microstructure ensemble n is larger than or equal to 200 for all computations: it is chosen, as given in Table 2, to achieve a small confidence interval for the average stress fields used for the parameter identification; see Sect. 4.3. We compute δ using (2), (10), and the area of the disc. In Fig. 18, the apparent material parameters are displayed for different area fractions A A and contrasts c using discs as grains. We also provide the RVE size as given in [37,42] (black circles connected by a black curve) when aiming at a relative error in the material parameters of 1 percent in the approach of Kanit et al., for a certain contrast c and using n = 200 microstructures samples.
We clearly see that a bias is introduced for small volume sizes, i.e., for small δ.
We also see that the asymptotic behavior of both parameters λ and μ is reached earlier (in terms of δ) for lower contrast; this is also reflected by the RVE sizes δ RVE (λ) and δ RVE (μ). For high contrast, i.e., c = 100, see Fig. 18 (right), we see clearly that the asymptotic behavior is reached faster for low area fraction (i.e., A A ≤ 0.4) as well as high area fraction (i.e., A A ≥ 0.8) than for an area fraction around 60 percent (i.e., 0.4 < A A < 0.8); this is, again, also reflected by the RVE size.
Especially for high contrast the change in the apparent material parameters is clearly visible. We can summarize that for small contrasts the RVE size tends to be smaller than for high contrasts and for very high and low area fraction the RVE tends to be smaller than for an area fraction around 50 percent.
From the results in Fig. 18, for simplicity, we choose, for discs, an RVE size of δ = 64, except for c = 100 and A A = 0.6 and A A = 0.7, where we use δ = 96, which is combined with a microstructure ensemble size of 200; see Table 2.

Effective parameters as a function of the area fraction -for different contrast
In Fig. 19 Table 2 In both diagrams, also the Voigt and Reuss bounds are shown, which are valid for linear elasticity. We see that our parameters are within these bounds.
In the bottom of Fig. 19, we show the same data as in the top plots, but here the effective parameters λ * and μ * have been normalized, i.e., the y-axis now has the range [0, 1]. Here, we see that for higher contrast c the deviation from the linear mixture rule, i.e., the Voigt bound, is significantly larger.
It is, indeed, well known in literature that, for different contrasts, the effective material parameters show a different dependency on the area fraction of the stiff inclusions and on the size of the RVE, cf. [46,70,74]. As a result, standard bounds (such as [16]) are not always useful, here.
In [86], the authors investigate the bulk modulus of a linear elastic material using a parametrized nonlinear power-law. They find that the effective bulk modulus decreases when the nonlinearity of their power-law increases, especially for larger area fractions of the inclusions phase.

Grain type and the size of the representative volume element
We also consider the grain type (and the specific boundary length L A ) and the size of the RVE. We use D, R 1 and R 3 as grain types. Note that we maintain the same grain resolution for all cases; see Fig. 20. In Fig. 21, the apparent material parameters are displayed for Boolean models with different area fractions and using discs and rectangles with different aspect ratio as grains.
Let us first note, that we observe that the macroscopic material behavior is stiffest for the elongated rectangles R 3 . We observe that the apparent material parameters seem to reach the asymptotic behavior earlier for rectangles R 1 than for discs and earlier for elongated rectangles R 3 than for  Fig. 17 Fit of the D 2 λ (δ) and D 2 μ (δ) over δ for discs and the apparent material parameters λ and μ using the approach by Kanit et al. [42]. For the fit only the blue data points were used. The influence of biased data on the fit can be derived from the distance between the red crosses and the curve of the fit for small δ. This is clearly visible in the plots for high contrast. The integral range δ 2 and the exponent α computed from this fit are shown in Table 4 Table 2. The RVE sizes δ RVE (λ) and δ RVE (μ) resulting from the approach by Kanit et al. [42] are also shown assuming the reference parameters rel = 0.01 and ensemble size n = 200   rectangles R 1 . It seems that the higher aspect ratio of the rectangles leads to earlier convergence of the parameters in terms of δ. The overlapping elongated rectangles indeed form connected clusters, stiffening the structure, earlier than, e.g., discs. We, again, note that the asymptotic behavior is reached slightly earlier for μ than for λ. It can also be observed that for discs and an area fraction around 60 percent (i.e., 0.4 < A A < 0.8) μ is larger than λ. For R 3 the opposite is the case.
With respect to the RVE sizes according to Kanit et al. [42], we see in Fig. 21 that δ = 21 is sufficient for the elongated rectanges R 3 ; see also Table 7. For the rectangles R 1 , using δ = 38.8, we do not fulfill the RVE criterion according to our reference choice of parameters (relative error rel = 0.01 and an ensemble size n = 200) for A A = 0.6. This is also reflected by the width of the confidence intervals for the stresses. Therefore, for this case we have used δ = 48.
From these results, we are confident to choose, except for the case with R 1 mentioned above and two cases with D (see Sect. 5.6), RVE sizes of δ = {64; 38.8; 23.4} for D, R 1 and R 3 , which are to be combined with microstructure ensemble sizes of 200; see Table 2. Note that these sizes δ arise from an image resolution of N 2 = 512 2 and a grain resolution of r N = 50 for the given grain types.

The shape of the grains
It is well known that the shape of the inclusions of the microstructures has an influence on the effective material properties of composite materials; cf. [12,44,68]. For each Boolean model, using different area fractions and using discs D and rectangles R 1 , R 3 , we construct 200 microstructures samples with δ = {64; 38.8; 23.4} for the given grain types.
We choose, according to Sect. 5.8, for discs D with area fractions A A = 0.6 and A A = 0.7 where we use δ = 96 and rectangles R 1 with area fraction A A = 0.6 where we use δ = 48. Again, we choose an ensemble size of n = 200. Making these exceptions, we ensure that we use large enough VEs; see Table 7 for the sizes of the RVEs.
In Fig. 2, we have shown examples of the microstructure samples used to compute the effective macroscopic material response. Note that the area of the discs and rectangles are chosen to be identical. However, the boundary length to area ratio L A /A A differs. In Fig. 22, we show the effective material parameters (top) and the normalized effective parameters  Fig. 21 The apparent material parameters of Boolean models using D, R 1 and R 3 as grain types for a contrast c = 100 are displayed. The microstructure ensemble sizes and the relative error for the apparent stresses are displayed in Table 2. We also show the RVE sizes δ RVE (λ) and δ RVE (μ) resulting from the approach by Kanit et al. [42] Table 3 The relative errors of rel (λ) and rel (μ), see (9), for c v = 1 are displayed for the computations in Figs. 24 and 25 (bottom) and observe that grain types with a low L A /A A ratio exhibit a softer macroscopic behavior.

Unstructured meshes
Next, we consider discs as grains where the diameter follows a log-normal distribution. Note that we do not have a maximum or minimum grain size in this case.
Since grains can now be large or small and, also, grains of different sizes displayed with structured meshes defeat the purpose of introducing a grain resolution, we move to unstructured meshes. To generate unstructured meshes for microstructure samples from the Boolean model, we use Mesh2D [11], see Fig. 15. Fig. 23 The apparent material parameters (bottom) are displayed for an increasing number of discs used for the different configurations as displayed in Fig. 12. The number of discs on the y-axis is given per dimension, i.e., #discs= 64 means that 4096 discs were placed on the unit square. Here, unstructured meshes were used (top). Identification with uniaxial-x/y tension test data  Fig. 24 The apparent material parameters of Boolean models using discs with log-normal radii with c v = {0; 0.125; 0.25; 0.5; 1} for a contrast c = 100 are displayed as a function of δ. Note, that the peak in the curves for A A = 0.2 and δ = 40 is an outlier, caused by a sin-gle disc in a microstructure sample that covers the entire domain. The microstructure ensemble sizes are determined heuristically as described in Table 2. We also show the RVE sizes δ RVE (λ) and δ RVE (μ) resulting from the approach by Kanit et al. [42] For the mesh generator, we use a minimum of 20 points to approximate the boundary of the grains. This is increased further for larger grains.

Sanity check for unstructured meshes
As a sanity check, we first consider regular patterns of discs as inclusions. We compare three configurations, i.e., centered, for the case of a single disc. By increasing the number of discs, we also increase δ, and we observe that the apparent parameters λ and μ converge at the same parameter value for large δ; see Fig. 23 (bottom).

Discs with log-normal radii
Finally, we investigate the influence of variability in the grain size. We investigate the effect of log-normal radii on the material parameters on the macroscale. We introduce variance in the distribution of the discs radii R ∼ LN (μ, σ ) using the coefficient of variation c v , where E(R) = e μ+σ 2 /2 is the mean and D R = E(R) e σ 2 − 1 is the standard deviation of R.
In Fig. 24, we choose volume element sizes δ = {8; . . . ; 96} and study the apparent material parameters for discs of constants size (c v = 0) and discs with c v = {0.125; 0.25; 0.5; 1} for constant area fraction A A . Note, that for c v = 0.5 and A A = 0.7 and A A = 0.8 we choose δ = 128 as well. We observe that for discs with high c v the asymptotic behavior is reached significantly later than for discs with constant size. We also have lower values for the material parameters for area fractions A A > 0.5.
Note, that the peak in the curves for λ and μ for A A = 0.2 and δ = 40 is caused by an outlier, i.e., a single disc in a microstructure sample that covers the entire domain, skewing the stress distribution of the microstructure ensemble, and resulting in a very large relative margin of error of the stress rel ( P 11 ) = 0.14.
In Fig. 25 the effective material parameters (top) and the normalized effective material parameters (bottom) are displayed for different c v . It holds that the macroscopic behavior is softer for high variance, i.e., up to 20 percent for c v = 1.
Note that by δ = 96 (see Fig. 25), we refer to the average since, in this context, δ is a distribution. As a result, for microstructure samples where very large discs occur, we may overestimate the stiffness as a result of the bias from the boundary conditions. On the other hand, very large discs lower the probability of connected clusters of inclusions (see Fig. 15), therefor softening the macroscopic material as we observe in Fig. 25.
With respect to the RVE size according to Kanit et al. [42], we see in Fig. 24 that using δ = 96 is not sufficient to fulfill the reference RVE criterion for our choice of parameters (relative error rel = 0.01 and an ensemble size n = 200) for two computations with c v = 0.5 and six computations with c v = 1. For the first cases we therefore use δ = 128 to fulfill that criterion.
For the cases with c v = 1, using δ = 96, we compute relative errors rel (λ) ≤ 0.03 and rel (μ) ≤ 0.025 in the material parameters, see Table 3, which may limit our interpretation of Fig. 25.
The relative margins of error of the stress are displayed for δ = 96 in Table 10. The results are consistent with the RVE size according to Kanit et al. [42]. We observe that a high variance in the distribution of the size of the discs results in a relative margin of error that is up to a magnitude higher, compared to discs of constant size. As a result, larger ensemble sizes are required to obtain a small statistical error for discs with log-normal radii.

Conclusions
We have combined the Boolean model, a well established model in stochastic geometry, with Neo-Hooke hyperelasticity to construct bi-phase microstructures from which we determined the effective material response in a systematic study.
The material parameters of the macroscopic Neo-Hooke model were computed by parameter identification using the ensemble averages of the mean stress and strain of many microstructure samples.
We have considered the confidence intervals for the stresses in choosing the ensemble sizes.
We have also successfully applied the approach of Kanit et al. [42] (see also [37]) to compute RVE sizes for the material parameters for different contrast, volume fraction, grain types, and variability in the grain size.
A large size characteristic δ of the microstucture samples was needed for high contrast c in the material parameters and grains with a low ratio of boundary to area of the grains. A high contrast in the material parameters results in a high deviation of the effective parameters from the linear mixture rule. Grains with a higher ratio of boundary to area induce a stiffer macroscopic material response.
Using log-normal disc radii larger RVE and larger microstructure ensembles are required to maintain low systematic and random errors. For higher variability in the distribution of the disc radii, we observed softer macroscopic response compared to discs of constant size. Funding Open Access funding enabled and organized by Projekt DEAL.
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://creativecomm ons.org/licenses/by/4.0/. Table 4 The fitted parameters δ 2 and α are displayed for different contrast and different types of grains and for both material parameters     Table 6 The effective parameters λ * and μ * are displayed for the computations in Figs. 19 and 22  Table 9 The fitted    Table 11 The effective parameters λ * and μ * are displayed for the computations in Fig. 25. Note, that this data was computed using unstructured meshes