Homogenization of Random Porous Materials With Low-Order Virtual Elements

A fast statistical homogenization procedure (FSHP) based on virtual element method (VEM)—previously developed by the authors has been successfully adopted for the homogenization of particulate random composites, via the definition of the representative volume element (RVE), and of the related equivalent elastic moduli. In particular, the adoption of virtual elements of degree one for modeling the inclusions provided reliable results for materials with low contrast, defined as the ratio between mechanical properties of inclusions and matrix. Porous media are then here described as bimaterial systems in which soft circular inclusions, with a very low value of material contrast, are randomly distributed in a continuous stiffer matrix. Several simulations have been performed by varying the level of porosity, highlighting the effectiveness of FSHP in conjunction with virtual elements of degree one.


Introduction
In recent times, the scientific community paid great attention to the influence of inherent uncertainties on behaviour of systems and recognized the importance of stochastic approaches to engineering problems [2]. Engineering experience has shown that uncertainties are involved not only in the assessment of loading but also in the material and geometric properties of engineering systems [3,4]. An example may be the mechanical characterization of complex realistic materials, such as heterogeneous media exhibiting non-periodic internal micro-structure. In this framework, statistical computational methods may be useful to evaluate their overall mechanical properties [5][6][7][8][9].
To this material typology belong porous solids: natural or artificial materials optimally designed for the load they carry and the environment in which they exist, such as porous rocks, cork, bones, sponges, aerogel, porous silicon and many others, that have great diffusion in engineering applications. They are, by their very nature, heterogeneous materials characterized The outline of the paper is as follows. In Section 2 the FSHP is detailed and the basic assumption of VEM are reported, in particular for the case of virtual elements of order 1. Section 3 provides the results of parametric analyses on random porous materials, modelled as two-phase material with soft inclusions. Attention is paid to the definition of the RVE size and to the sensitivity to porosity: how RVE size and mechanical properties change at varying of porosity. In Section 4 some final remarks are presented.

Fast Statistical Homogenization Procedure with VEM
In this Section, the steps of the automated so-called Fast Statistical Homogenization Procedure (FSHP) [1], based on the statistical homogenization procedure previously developed in [23] and briefly described in Section 1, are detailed. FSHP allows to compute overall homogenized mechanical parameters, in particular the anisotropic elasticity tensor. Boundary Values Problems (BVPs) are solved at microscopic level, where each material phase is characterized by linear elastic isotropic behaviour, applying both Dirichlet and Neumann boundary conditions. At the macroscopic level, a general anisotropic 2D stress-strain relation is derived applying the Hills macrohomogeneity condition. Afterwards, the virtual element of degree one adopted in this work is described.

Constitutive model
We consider a two-dimensional framework and describe, at the microscopic scale, the heterogeneous porous material as a bi-phase system. In two-phase materials it is useful to define the material contrast as the ratio between the elastic moduli of inclusions, E i , and matrix, E m , c = E i /E m . When 0 < c < 1 (c = 1 being the case of a homogeneous material) inclusions are softer than the matrix and we refer to low contrast materials, that are suitable to properly represent porous media [11]. The pores are modelled as circular inclusions of diameter d dispersed into continuous matrix ( Fig. 1(a)). Furthermore a scale parameter δ = L/d is introduced, that is the ratio between the side of a square test window L, and the diameter d of the inclusions. The analysis are developed within a kinematic linearised framework.
In order to perform homogenization, we consider two levels of the description, the microscopic and the macroscopic level, with governing equations that are formally the same except than for the constitutive law that is not 'a priori' defined at the macroscopic level, but directly descends from the lower level as result of the homogenization procedure. In the following, lower case letters always refer to the microscale, while upper case letters to the macroscale. We assume that at the lower level each material phase is characterized by linear elastic isotropic behaviour with the stress-strain relations: where ε i j and σ i j (i, j, k = 1, 2) are the components of the microstrain and microstress tensors, λ and µ are the Lamé constants, and δ i j (i, j = 1, 2) the Kronecker symbol. At the macroscopic level, the general anisotropic 2D stress-strain relations, read: where E i j , Σ i j (i, j = 1, 2) are the components of the macrostrain and macrostress tensors and C i jhk (i, j, h, k = 1, 2) the components of the macroscopic elastic tensor, corresponding to the equivalent moduli obtained via a homogenization procedure based on the Hill's macrohomogeneity condition [38]: which states the equivalence between the stress density power in the macroscopic model and the mean power evaluated at the microscopic level within a domain (test window), B δ , occupying a region of area A δ .
The set of Dirichlet and Neumann boundary conditions (BCs) to apply, at the microscale, on the boundary ∂B δ , for the solution of the BVPs, directly derive from the fulfilment of the macrohomogeneity condition (3). The Dirichlet BCs are: u i being the components of the displacement vector and x j (i, j = 1, 2) the j-th coordinate of the generic point on the boundary, ∂B δ , with respect to a reference system with origin in the geometric center of the test window B δ ; while the Neumann BCs are: where t i are the components of the traction vector and n j the component of the outward normal to ∂B δ (i, j = 1, 2).

Statistical Homogenization
The homogenization procedure for defining the constitutive response of random heterogeneous material requires the definition of the size of a Representative Volume Element (RVE) larger than the microscale characteristic length, corresponding to the diameter of inclusions, d.
According to the approach presented in [23], based on the approach proposed in [22,25], the presented procedure requires the statistical definition of a number of realizations, called Statistical Volume Elements (SVEs), of the possible microstructure, sampled in a Monte Carlo sense, which allows determining series of scale-dependent upper and lower bounds for the overall elastic moduli and to approach the RVE size, δ RV E , using a statistical criterion to stop the procedure when the results in terms of average elastic moduli do not change within the selected tolerance interval.
All the steps of the homogenization procedure are completely integrated in the FSHP. The entire computational stochastic homogenization procedure is structured as follows and resumed in the Flow chart in Fig. 2.  Step 1 Input: set the nominal volume fraction of the medium (ρ ≤ 40%) as the ratio between the total area of the inclusions and the area of a test window. The particles are disk-shaped inclusions, with diameter d, randomly distributed on the base of a Poisson point field process with 'hard-core' distance, s, between the particles, fixed in order to avoid very narrow necks between inclusions. Define the dimensionless scale factor δ = L/d. Fix the mechanical parameters of each phase: Young modulus of inclusions and matrix, E i and E m ; Poisson coefficient, ν i and ν m . Set the minimum number of simulations for convergence, N lim , and a tolerance parameter, Tol, based on data dispersion, as defined in above.
Step 2 Input: initialize the window size, L = L 0 , and number of simulations, N = N 0 .
Step 3 Realizations: for each window size determine the number of inclusions (Poisson random variable) via simulations by means of Knuth's algorithm, where the Poisson's distribution is characterized by the parameter p δ = (ρL 2 )/[π(d/2 + s/2) 2 ]. Then simulate random dispositions of disks' centres, according to the uniform distribution, that are the realizations (SVEs) of portions of the random medium. Each realization is supposed to be independent from any previous one. An example of realization of the miscrostructure and the related generated mesh is shown in Fig. 1(a).
Step 4 Generate/Solve: for each SVE, generate the relative mesh ( Fig. 1(b)) and solve both the Dirichlet (Equation 4) and Neumann (Equations 5) BVP, and compute the homogenized constitutive parameters.
Then compute N i = (1.96 σ(K δ )/(K δ Tol) 2 , which ensures that the confidence interval of the average homogenized constitutive parameter set at 95%, evaluated over the normal standard distribution, is within the tolerance allowed, Tol.
Repeat Steps 3-4 until N i < N lim .
Step 6 Checking: if the number of realizations necessary for ensuring the requirement at Step 5 is small enough, stop the procedure. We choose as the number of realizations necessary the most unfavourable number between those obtained by solving BVPs of Neumann or Dirichlet. Otherwise choose an increased value of δ and go to Step 3.
The fulfilment of the requirement at Step 6 means that the values of the homogenized constitutive coefficients are distributed around their averages with a vanishing variation coefficient, and that the RVE size is achieved. The effective constitutive moduli are finally estimated as the mean values between the Dirichlet (upper) and Neumann (lower) bounds at the convergence window (RVE).
The statistical convergence criterion adopted is based on the 95% over the Normal Standard distribution, which provides the number N at which is possible to stop the simulations for a given window size δ. When this number is small enough, the average values of the effective moduli converge and the RVE size is achieved. This circumstance also corresponds to reaching the minimum window size δ RV E for which the estimated homogenized moduli remain constant, within a tolerance interval less than 0.5% for both the Dirichlet and Neumann solutions. The minimum number of simulations, N lim , and the tolerance parameter, Tol, are chosen in order to define a narrow confidence interval for the average and to obtain a reliable convergence criterion. The choice is discretionary, values are assumed depending on the data dispersion. All these details are specified in [23].
The FSHP permits the automatic cutting of the inclusions over the limit of the windows' edges, accounting for the presence of inclusions that randomly cross the windows' edges as required by the randomness of the medium. It is worth noting that to neglect the presence of inclusions that intersect the windows' edges, a less realistic hypothesis widely used in literature, provides results significantly different from the results obtained taking into account cut inclusions at the windows' edges [36]. Furhermore, the mesh corresponding to each realization of the microstructure has been optimized using VEM methodology, that permits to adopt single virtual element for the inclusions (reduction of degrees of freedom) and triangular virtual elements for the matrix ( Fig. 1(b)). In order to take into account the stress gradients in the so-called 'hard core regions', the mesh is fine near the inclusions and coarse away from the inclusions. The FSHP allow us to very efficiently solve the series (hundred) of BVPs and to rapidly converge to the RVE solution.

Virtual Element Method
In this section we recall the weak formulation of the classical linear elastic problem and describe the related virtual element space of degree 1 [42,43] as well as the construction of the bilinear form resulting from the weak form.
In the following we use the Voigt notation, suited to the proposed formulation.

The linear elastic problem
Let Ω be a two-dimensional domain, and let Γ be its boundary. We consider the deformation problem of a linearly elastic body subjected to the volume force, represented by the vector f f f ∈ L 2 (Ω) 2 where ∂ (·) denotes the partial derivative with respect to the (·)-coordinate.
The weak form of the linear elastostatic problem, provided by the virtual work principle, reads: where: where C C C is the 3 × 3 elastic tensor.

Virtual element formulation
In this subsection we introduce the virtual element discretization used in the following simulations. In order to approximate the solution of the problem (7) we consider a decomposition T h of the domain Ω into non overlapping polygonal elements E. In the following, we denote by e the straight edges of the mesh T h and, for all e ∈ ∂E, n n n i denotes the outward unit normal vector to e i ( Fig. 3(a)). The symbol n e represents the number of the edges of the polygon E. Let k be an integer ≥ 1. Let us denote by P k (Ω) the space of polynomials, living on the set Ω ⊆ R 2 , of degree less than or equal to k.
By the discretization introduced, it is possible to write the bilinear form (7), as in the finite element methodology, in the following way: The discrete virtual element space, V V V h , is: where the local space, V V V h|E := V h|E 2 . For the virtual element of degree k = 1 V h|E is defined as: The dimension of the space V V V h|E then is: We can observe that, in contrast to the standard finite element approach, the local space V V V h|E is not fully explicit. Moreover v v v h is a polynomial of degree 1 on each edge e of E and globally continuous on ∂E.
The problem (7) restricted to the discrete space V V V h becomes: where a h (·, ·) : V V V h ×V V V h → R is the discrete bilinear form approximating the continuous form a(·, ·) and, < f f f , v v v h > is the term approximating the virtual work of external load. In this work the body force is null and then we do not report any detail for the implementation. The discrete bilinear form is constructed element by element as By the above definition of the local space (11) the following important observations can be made: The related degrees of freedom for the space V V V h|E are 2n e point-wise values v v v h (υ i ) i = 1, 2, . . . , n e , where υ i is the i-th corner of E.

Projection operator and construction of the stiffness matrix
In accordance to [51] we define the projector operator of the strain as: In particular the projector operator respects the following property of orthogonality: As in the finite element method, we define the vectors The matrices N N N u and N N N V are where φ i is the standard i−th shape function. The projection operator in matrix is: with Π ∈ R 3×m . By putting Equations (17)(18)(19) into the Equation (16) we obtain: Integrating by part the right hand side we obtain: where N N N E contains the components of the outward normal n n n = {n 1 , n 2 } T N N N E = n 1 0 n 2 0 n 2 n 1 .
After some algebra: Finally the local projector operator Π is: where: We can observe that the matrix G G G ∈ R 3×3 is computed by knowing the area of the elements | E |: and the matrix B B B ∈ R 3×m using the degrees of freedom of the element: (n 1 1 | e 1 | +n n e 1 | e n e |)/2 0 · · · (n n e 1 | e n e | +n (n e −1) 1 | e (n e −1) |)/2 0 0 (n 1 2 | e 1 | +n n e 2 | e n e |)/2 · · · 0 (n n e 2 | e n e | +n (n e −1) 2 | e (n e −1) |)/2 (n 1 2 | e 1 | +n n e 2 | e n e |)/2 (n 1 1 | e 1 | +n n e 1 | e n e |)/2 · · · (n n e 2 | e n e | +n where | e i | is the length of the i−th edge of the element E and n i j (i = 1, · · · , n e ; j = 1, 2) is the i-th outward normal of components j.
In the virtual element strategy we approximate the bilinear form of Equation (14) as follow where the first term of the right hand side is the consistent term and the second term is the stabilization term. For computing the consistent term M M M we use Equation (19) and Equation (8) where If the elastic tensor is constant, it is possible to rearrange Equation (30) and obtain: Regarding the stabilization term, we adopt the choice introduced in [42,43]. First of all, we introduce the scaled coordinates (ξ, η) as: where (x E , y E ) are the coordinates of the centre of E and h E its diameter.
The stabilization term are define as follow where α ∈ R and in this work is set equal to 1/2. The matrix D D D contains the scaled coordinates at the nodes of the element: 3 Numerical analyses In this section, the Fast Statistical Homogenization Procedure (FSHP) based on the low order virtual elements is applied to the analysis of porous materials, modelled as a two-dimensional two-phase material in which circular soft inclusions are randomly embedded in a stiffer continuous matrix.
FSHP with virtual element of k = 1 is particularly suitable for analysing low contrast materials [1], and in this paper a very low value of contrast (c → 0) is adopted, with the purpose of simulating a material with randomly distributed voids. Referring to the properties adopted in [11], an high value of Poisson coefficient has been adopted for the inclusions. All the mechanical properties are reported in Table 1, where E m and ν m are Young and Poisson modulus for the matrix, and E i and ν i are Young and Poisson modulus for the inclusions, respectively.  Table 1. Mechanical properties adopted The procedure is applied to the definition of the RVE size in case of random porous materials and to the evaluation of their mechanical properties. Initially the procedure is used keeping constant the porosity, defined as the volume fraction of pores, at a value ρ = 40%, for describing how the RVE size, δ RV E , is identified. Subsequently, exploiting the opportunities provided by FSHP, several parametric analyses are performed by varying the porosity, ρ, in the range 0% ÷ 40%. The purpose is twofold: on one hand for identifying δ RV E and its changes in relation of porosity; on the other hand, to evaluate the sensitivity to porosity of the mechanical properties of homogeneous equivalent continuum materials.

RVE definition
Homogenization of random porous materials, modelled as a very low contrast material, requires an high number of realizations and the adoption of large windows to detect the RVE size, δ RV E , and determine the overall elastic moduli [37]. FSHP consists in solving the Boundary Value Problem (BVP) on a huge number -hundreds of analyses -of test windows of increasing size, δ. BVPs are solved by applying both Dirichlet and Neumann BCs, providing upper and lower bounds for the equivalent elastic moduli respectively. As the windows size increases, the number of simulations needed to obtain the convergence of the results decreases. By referring to Section 2.2 we recall that convergence of the results is reached when the average value of an estimated mechanical parameter does not significantly changes after N simulations. The homogenized elasticity tensor is anisotropic by definition, however, the pores being of circular shape and randomly distributed, after several simulations the homogenized elastic tensor moduli has been found to be substantially isotropic. Then in the following we select the bulk modulus, K = 2C 1122 +C 1212 as representative material parameter to identify.
In Fig. 4, the average value of the bulk modulus obtained after N simulations K n is normalized versus the value obtained after N − 1 simulations, both for Dirichlet BCs (Fig. 4(a)) and for Neumann BCs (Fig. 4(b)). Two different windows of increasing size δ are considered as example: δ = 20, in blue solid line, and δ = 30, in red solid line. It is possible to notice that the value of the ratio K n /K n−1 undergoes less variations as the number of simulations N increases, and that, as expected, the number of simulations needed to obtain the convergence decreases as the windows size δ increases. In the reported example, the convergence is reached for N 55 for δ = 20 and N 25 for δ = 30 in Dirichlet BVP, while N 125 for δ = 20 and N 45 for δ = 30 are necessary in the case of Neumann BVP. After reaching the convergence for a given δ, FSHP goes on and larger windows' size, δ, are analysed until the number of simulations needed to obtain convergence, N, is small enough (Step 6 in 2.2) Based on the criterion of the minimum number of simulations needed, the RVE can be reached also for material with very low contrast, such as the porous materials here considered. Note that δ RV E is then identified as the smallest window size needed to reach the convergence for the BVPs, considering the most unfavourable case between Dirichlet and Neumann BCs. The results obtained by applying Dirichlet and Neumann BCs provide the same trend, however it is possible to notice that the number of simulations needed is smaller in Dirichlet BVP solution. The difference in term of convergence between the two BVPs has been highlighted in [5,22]. In particular, it has been showed that in the case of material with high contrast, where inclusions are stiffer than matrix, Neumann BVP reaches the convergence faster than Dirichlet BVP, while the opposite occurs for low contrast materials [1,37]. In the case of porous media, the most unfavourable case is Neumann BVP, meaning that larger windows are needed to reach the convergence value with respect to Dirichlet BVP. In Fig. 5, the number of simulations needed to obtain the convergence for increasing value of windows size δ is reported. Red solid lines are referred to the solution of Dirichlet BVP, while blue solid lines to the solution of Neuman BVP; the same colors are used for Dirichlet and Neumann in the following figures. Dirichlet BCs provides convergence results with a low number of simulations, starting from windows of size δ ≥ 25, while larger windows, δ > 35, are needed when adopting Neumann BCs. The identification of δ RV E always requires the solution of the Neumann BVP, resulting, as mentioned above, the most unfavourable case.  6 reports the average value of the bulk modulus, K, normalized with respect to the value estimated at the RVE, K RV E , versus the number of simulations performed for different window sizes, δ, providing a picture of the statistical convergence criterion adopted. K RV E is evaluated as the mean value between the two values obtained by solving Dirichlet and Neumann BVPs at the convergence window δ RV E . δ RV E is identified for a window size δ 37, obtained after about N 20 simulations (Fig. 5). As expected, widening the windows size provides decreasing K in Dirichlet BVP solution ( Fig. 6(a)) and raising K in Neumann BVP solution ( Fig. 6(b)). Concerning the convergence trend, as δ increases the results of Dirichlet and Neuman BVPs tend to get closer to each other, as shown in Fig. 7(a), where the average value of the homogenized bulk modulus normalized with respect to the value at δ RV E , K/K RV E , is plotted versus the scale parameter δ. Even if there is a mismatch between the two solutions, for values of δ ≥ 37 they become parallel and δ RV E is identified.
The differences in terms of convergence trend between Dirichlet and Neumann BVPs solutions also depend on the different dispersion of the results, as shown in Fig. 7(b), where the Coefficient of Variation, CV , is plotted versus δ. Note that the convergence is obtained for values of CV ≤ 0.05 and this occurs for δ ≥ 25 in Dirichlet BVP and for δ ≥ 37 in Neumann BVP.  The results of the parametric analyses performed are reported in the following figures. As above, red color always refers to the solution of Dirichlet BVP, while the blue one to the solution of Neumann BVP. The graphs in Fig. 9 provide a synthetic response of the mechanical parameters variation for different values of porosity. Fig. 9(a) summarizes the results of the parametric analyses in term of bulk modulus K versus the porosity ρ. As previously obtained, for increasing δ, the two bounds tend to become closer to each other. The value of K decreases in Dirichlet solutions while raises in Neumann solutions, as shown in the zoom area of Fig. 9(a) for ρ = 25, where the values of k for increasing δ are plotted both for Dirichlet (circles) and Neumann (squares) BVPs. As underlined in Section 3.1, Dirichlet BVP reaches convergence for δ smaller while Neumann BVP requires larger δ and results are more dispersed.
For each value of ρ, the lowest value of K provided by Dirichlet BVPs and the higher value obtained by Neumann BVPs are achieved at δ RV E . The area included between Dirichlet and Neumann BVPs, marked by grey lines pattern, represents the envelope between the two bounds: K RV E for increasing ρ is obtained as the mean value between the two bounds. As expected, the value of K decreases as the porosity increases [48]. Apparently the two bounds tend to become closer as ρ increases, however the percentage difference between the bounds and the mean value remains almost constant. Fig. 9(b) provides the results of the parametric analyses in term of Poisson ratioν versus porosity ρ. Also in this case, ν decreases as the porosity increases, as reported in scientific literature [49,50]. It is interesting to point out that in this case, differently from what generally occurs, Neumann BVP solutions provide higher values of ν with respect to Dirichlet BVP, becoming the upper bound of the solution. This is probably due to the fact that Neumann BCs are weaker then Dirichlet ones, where transversal contraction is avoided, so that transversal strain is bigger and consequently Poisson's coefficient too.
(a) Bulk modulus K for different density ρ (b) Poisson modulus ν for different density ρ It is interesting to look out at the RVE size obtained for different values of porosity, (Fig. 10). In the case of ρ = 0 the material is homogeneous, starting from ρ ≥ 0 it is necessary to determine a RVE and, as expected, as porosity increases bigger δ RV E are needed. However, it has been found that the δ RV E rapidly increases in the range of porosity 1 ÷ ρ ÷ 20, while slowly decreases up to ρ = 25 and then remain almost constant in the range 25 ÷ ρ ÷ 40. The results in terms of RVE size, δ RV E , are synthesized and reported in Fig. 10, where the values of δ = δ RV E are plotted versus the different values of porosity ρ considered. As previously noticed, Dirichlet BVP need smaller windows with respect to Neumann BVP, that has to be considered to define δ RV E for porous material.
The differences in terms of convergence trend between the two BCs depend on the different dispersion of results, as shown in Fig. 11, where the Coefficient of Variation CV is plotted versus several values of porosity ρ for Dirichlet ( Fig.  11(a)) and Neumann BVP (Fig. 11(b)).

Conclusions
In this work, a fast statistical homogenization procedure (FSHP) has been adopted for the homogenization of random porous materials. Porous material has been modelled as bi-phase material with very low contrasts: disk shaped soft inclusions randomly distributed within a base stiffer matrix. FSHP based on virtual elements of degree 1 allow us to perfectly model porous material [1] and permits to solve high number of simulations as required in homogenization techniques applied to random materials [23,36,37,40]. Furthermore, FSHP permits to analyse a series of materials with different porosity and to rapidly identify the relative RVE size and the overall effective moduli of the equivalent homogeneous material.
Some final remarks may be drawn: -FSHP provides reliable results for the homogenization of porous materials.
-VEM of degree 1 perfectly fits the behaviour of porous materials: stress and strain are constant in low order virtual elements, however this does not have influence on homogenization results in case low contrast materials, as shown in [1]. -The Representative Volume Element (RVE) as a function of the porosity has strongly non-linear behaviour. Next step of this research is to improve the accuracy of the results by adopting virtual elements higher degree for studying material with hard inclusions dispersed in softer matrix. Moreover, the proposed method will be implemented in order to better describe the behaviour of micro-structured material made of particles of different shape and orientation by adopting enriched non-classical continua descriptions, such as Couple-Stress and Cosserat, [40].