STATISTICAL HOMOGENIZATION OF RANDOM POROUS MEDIA

In recent times, the scientific community paid great attention to the influence of inherent uncertainties on system behavior and recognize the importance of stochastic and statistical approaches to engineering problems [21]. In particular, statistical computational methods may be useful to the constitutive characterization of complex materials, such as composite materials characterized by non-periodic internal micro-structure. Random porous media exhibit a microstructure made of randomly distributed pores embedded into a continuous matrix. They can be modelled as a bi-material system in which circular soft inclusions (pores) with random distribution and variable diameters are dispersed in a stiffer matrix. A key aspect, recently investigated by many researchers, is the evaluation of appropriate mechanical properties to be adopted for the study of their behaviour. Differently from classical homogenization approaches, in the case of materials with random microstructure it is not possible to ‘a-priori’ define a Representative Volume Element (RVE), this being an unknown of the problem. Statistical homogenization procedures may be adopted for the definition of equivalent moduli able to take into account at the macroscale the material properties emerging from the internal microstructure with random distribution [26]. Here, a Fast Statistical Homogenization Procedure (FSHP) based on Virtual Element Method (VEM) approach for the numerical solution – previously developed by some of the authors [13] has been adopted for the definition of the Representative Volume Element (RVE) and of the related equivalent elastic moduli of random porous media with different volume fraction, defined as the ratio between mechanical properties of inclusions and matrix. In particular, FSHP with virtual Elements of degree 1 [2] for modelling the inclusions provides reliable results for materials with low contrast.


INTRODUCTION
The characterization of porous materials is increasingly attracting the attention of researchers and engineers due to the widespread use of this class of materials in different fields and industries, including aerospace, automotive, energy, construction, electronics and biomedical. Focusing the attention on porous metals, porous ceramics and polymer foams, the main common features are low density, large specific surface and a range of novel properties ranging from physical, mechanical, thermal, electrical up to acoustical fields [10]. In addition, a distinctive characteristic of such heterogeneous materials is the random distribution of voids within a dense solid. This results in a composite two-phase material whose overall elastic properties depend both on the geometrical nature of the pores -shape and size -and on the value of porosity. Different approaches have been proposed in literature to handle the study of porous materials, among others we refer to generalized continua in which voids or microcracks are modelled as additional degrees of freedom [22,28,27,20], multifield [14,6,7,8,19] and/or multiscale descriptions [25,26]. This paper is devoted to the investigation of porous materials within the framework of linear elasticity. A fast statistical homogenization procedure (FSHP) is adopted to grasp the global elastic behaviour of such a random composite, as in [13]. In FSHP the statistical procedure, proposed in [26], is automatized and integrated in a completely in house specifically developed code implemented to quickly and efficiently perform a high number of parametric analyses. The material is modelled considering disk shaped soft inclusions randomly distributed within a base stiffer matrix. A first order computational homogenization procedure is exploited within the very recent Virtual Element Method (VEM) [4,5,13,12] that allows us to reliably model porous material and to efficiently solve high number of simulations as required in homogenization techniques applied to random materials [15,1,3]. The main point of the procedure is to approach the so-called Representative Volume Element (RVE) using finite-size scaling of Statistical Volume Elements (SVEs). To this end properly defined Dirichlet and Neumann-type boundary value problems are numerically solved on the SVEs defining hierarchies of constitutive bounds. A set of materials with different porosity is analysed to characterize overall mechanical parameters of porous materials and to investigate their sensitivity to porosity. It emerges that, on the one hand FSHP provides reliable results for the homogenization of porous materials. On the other hand, the choice of virtual elements of degree 1 is perfectly suitable to the case at hand, as shown in [13]. Moreover, both homogenized values of bulk modulus and Poisson coefficient decrease as the porosity increases.

FAST STATISTICAL HOMOGENIZATION PROCEDURE
We consider a two-dimensional linear elastic 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 [18]. The pores are modelled as circular inclusions of diameter d dispersed into continuous matrix. 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 homogenization procedure for defining the constitutive response of random heteroge-neous 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 [26], that is based on the approach proposed in [9,11], 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 so-called Fast Statistical Homogenization Procedure (FSHP), based on the statistical homogenization procedure previously developed in [26]. See [13] for details on FSHP procedure. The convergence criterion is fulfilled when the values of the homogenized constitutive coefficients are distributed around their averages with a vanishing variation coefficient. This means that the RVE size is achieved. The effective constitutive moduli are consistently estimated as the mean values between the Dirichlet (upper) and Neumann (lower) bounds at the convergence window (RVE). 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, T ol, 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 [26]. It is worth noting that 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 by taking into account cut inclusions at the windows' edges [25]. Furthermore, 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. 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 briefly recall the weak formulation of the classical 2D linear elastic problem and describe the related virtual element space [4,5], as well as the construction of the bilinear form resulting from the weak form. The vectorial notation is adopted in the following, that is suited to the proposed formulation. We consider a body immersed in the two-dimensional space R 2 , where the Cartesian coordinate system (O, x, y) is introduced. The body is subjected to the volume force, represented by the vector f ∈ (L 2 (Ω)) 2 , f = {f 1 , f 2 } T (within the standard Lebesgue space), and given boundary conditions. In the sake of simplicity we use homogeneous Dirichlet boundary conditions and consider the Sobolev space, V := (H 1 0 (Ω)) 2 , of the admissible displacement fields, represented by the vector v.
Furthermore we represent the strain, under the hypothesis of small deformations, as the vector ε = {ε 11 , ε 22 , ε 12 } T associated to the displacement field vector u = {u 1 , u 2 } T : where ∂ (·) denotes the partial derivative with respect to the (·)-coordinate.
The weak form of the linear elastic problem reads: where: and C = C(x) is the plane elastic tensor (uniformly positive) and possibly depending on the position vector x = (x, y) T ∈ Ω.
In order to approximate the solution of the problem (2) 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 i denotes the outward unit normal vector to e i ( Fig. 1(a)). The symbol n e represents the number of the edges of the polygon E, that coincides with the number of the element vertices. 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 (2), as in the finite element methodology, in the following way: The discrete virtual element space, V h , is: where V h|E := V h|E 2 and the local space V h|E is defined as By the definition of the local space (6), we can observe that, in contrast to the standard finite element approach, the local space V h|E is not fully explicit, in fact V h|E contain all the polynomials of degree ≤ k, plus other functions that, in general, will not be polynomials. Moreover v h is a polynomial of degree k on each edge e of E and globally continuous on ∂E.
The related degrees of freedom for the space V h|E (Figs. 1(b),1(c)) are: • 2n e point-wise values v h (υ i ) i = 1, 2, . . . , n e , where υ i is the i-th corner of E; • 2n e (k − 1) point-wise values v h (y e j ), where y e j , i = 1, . . . , k − 1 are the edge internal nodes; • k(k − 1) scalar moments of the unknown field over the element (i.e. 1 |E| v h , where | E is the area of the element), not associated with a specific location over E.

The global dimension of the space
The problem (2) restricted to the discrete space V h becomes: where a h (·, ·) : V h × V h → R is the discrete bilinear form approximating the continuous form a(·, ·) and, < f , v h > is the term approximating the virtual work of external load. The discrete bilinear form is constructed element by element as: The local stiffness matrix can be derived, consistently with [5], after introducing the local projector operator Π ∇ E : V h (E) → (P k (E)) 2 . The details on the construction of the local stiffness matrix are presented in [13], while are here omitted for the sake of brevity.

NUMERICAL SIMULATION: SENSITIVITY TO POROSITY
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 degree k = 1 is particularly suitable for analysing low contrast materials [13], and in this work 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 [18], 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.
By exploiting the opportunities provided by FSHP, several parametric analyses have been performed by varying the porosity, ρ, of random porous media in the range 0% ÷ 40% (Fig.  2). The purpose is twofold: on one hand for identifying the RVE size,δ RV E , and its changes in 10 −4 0.30 0.49 Table 1: Mechanical properties adopted Figure 2: Realizations of micro-structure for different levels of porosity relation of porosity; on the other hand, to evaluate the sensitivity to porosity of the mechanical properties of homogeneous equivalent continuum materials. The performed analyses allowed us to identify RVE size, δ RV E , and the corresponding material properties. As the homogenized material has found to be essentially isotropic, we focus the attention only on the bulk modulus, as representative material parameter. The results of the parametric analyses performed are reported in the following figures.
Basing on the works [16] and [23] we define, the following constitutive scaling measures are: where K D δ and K N δ are the average of the bulk moduli obtained by solving the boundary value problems (BVPs) by adopting Dirichlet and Neumann-type boundary conditions (BCs), respectively. Fig.(3(a)), provides a qualitative and quantitative information about the convergence trend of the solution,f K δ , by varying the window size, δ, and the pore density, ρ. The differences in terms of convergence trend between the two BC solutions depend on the different dispersion of results, as shown in Fig. 3(b), where the Coefficient of Variation, CV , is also plotted for several values of porosity ρ for Dirichlet and Neumann BVPs. Fig. 4(a) summarizes the results of the parametric analyses in term of convergence bulk modulus K, both for Dirichlet and Neumann BVPs, versus the porosity ρ.
As expected, the value of K decreases as the porosity increases [15]. 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. 4(b) reports the values of RVE size obtained for δ = δ RV E plotted versus the different values of porosity ρ. It is worth noting that 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 larger δ 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 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.

CONCLUSIONS
The present work is an application of the Fast Statistical Homogenization Procedure (FSHP) [13] to the homogenization of random porous materials. The model adopted within the FSHP framework is the bi-phase material in which disk shaped soft inclusions are randomly distributed in a stiffer matrix. FSHP uses the Virtual Element technique to model the inclusions with one element, this permit to reduce the number of degrees of freedom with consequently increasing the computational efficiency. FSHP permits to solve high number of simulations as required in homogenization techniques applied to random materials [25,26,24,17]. 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. Virtual Elements of degree one, adopted in this work, well fits the behaviour of porous materials. The approximation of constant stress and strain, using lower virtual element, in fact does not influence the homogenization procedure [13,12]. For porous materials the Representative Volume Element (RVE) as a function of the porosity has strongly non-linear behaviour and the RVE increases when the porosity increases, but for higher level of porosity the RVE slightly decreases and then remains constant up to the value ρ = 40%, which is the maximum value of porosity. Furthermore, as expected, the homogenized values of the elastic moduli decrease as the porosity increases.