Statistical Homogenization of Polycrystal Composite Materials with Thin Interfaces using Virtual Element Method

Polycristalline materials with inter-granular phases are modern composite materials extremely relevant for a wide range of applications, including aerospace, defence and automotive engineering. Their complex microstructure is often characterized by stochastically disordered distributions, having a direct impact on the overall mechanical behaviour. In this context, within the framework of homogenization theories, we adopt a Fast Statistical Homogenization Procedure (FSHP), already developed in Pingaro et al. (2019), to reliably grasp the constitutive relations of equivalent homogeneous continua accounting for the presence of random internal structures. The approach, combined with the Virtual Element Method (VEM) used as a valuable tool to keep computational costs down, is here successfully extended to account for the peculiar microstructure of composites with polycrystals interconnected by thin interfaces. Numerical examples of cermet-like linear elastic composites complement the paper.


Introduction
Modern composites materials are conceived to fulfil specific design requirements calling for the need of materials with increasingly high performances in terms of stiffness, strength, fracture toughness and lightness, among others. Special attention has been devoted in the last decades to advanced ceramics, namely composites with ceramic (CMC) or metal matrix (MMC) [1,2] which are gradually more used for a wide range of challenging applications, ranging from bioengineering, with the production of biocompatible ceramics for prosthesis and artificial organs [3], aerospace, defence, automotive [4] up to mechanical engineering for the production of wearing parts, seals, low weight components 10 and fuel cells, and cutting tools with distinguished thermo-mechanical and wear properties [5]. Often ceramic composites exhibit a peculiar microstructure, characterized by polycrystals interconnected by a second-phase interphase of small thickness in comparison to the grain diameter. Consistently with [6], different combinations of materials lead to four kinds polycrystalline materials with inter-15 granular phases: i) ceramic grains with brittle interphases; ii) ceramic grains and a metallic interphase; iii) metallic grains surrounded by ceramic interphases; iv) metallic grains surrounded by a layer of soft ductile material. In all cases, such materials show complex phenonema occurring at different scales of interest, ranging from the microscopic up to the macroscopic one, so that the elastic 20 properties of the homogenized composite material is strongly influenced by the random distribution of grains, by the role played by the intergranular layer, as well as by the value of the contrast, that is defined as the ratio between the elastic moduli of the inclusions and the matrix [7,8,9]. 25 The reliable evaluation of homogenized constitutive properties to be adopted for the macroscopic investigation of their behaviour pertains to the well-established homogenization theories, and is a long-standing problem of great interest for many researchers [10, 11,12,13,14,15,16], also with the purpose of designing innovative smart materials [17,18]. In this paper, we adopt a statistical ho- 30 mogenization approach [11,19] to derive the in-plane linear elastic constitutive properties polycristalline materials with thin layers interphases and to evaluate how they change in relation to contrast.
In the case of materials with random micro-structure, the lack of periodicity in 35 the microscopic arrangement makes the homogenization process more challenging, asking for special attention in the detection of the Representative Volume Element (RVE) along with the evaluation of homogenized constitutive moduli.
Several approaches have been proposed [20,21,22,23,24,25,26,27,14,13,28,29], also referred to non-classical continua [30,19,9]. 40 Among these models, we focus our attention on the possibility of approaching the RVE using finite-size scaling of intermediate control volume elements, named Statistical Volume Elements (SVEs), and proceed to homogenization (e.g. [31]). In this work we refer to [19], where a homogenization procedure consistent with the Hill-Mandel condition [32] has been coupled with a statisti- 45 cal approach, by which scale-dependent bounds on classical moduli are obtained using Dirichlet and Neumann boundary conditions (BCs) for solving boundary value problems (BVPs). This statistically based homogenization procedure has provided significant results [30,33,34], with particular reference to the debated problem of the convergence in the presence of materials with very low (or very 50 high) contrast [9,31,25,35].
Here, the procedure is to specifically account for the microstructural topology exhibited by polycrystalline composites with thin interphases, making the most of VEM to improve the computational efficiency. The result is a powerful numerical tool suitable for analyses of large portions of random microstrure, very important in the case materials with very high or low contrast [37]. In particular the strategy proposed is to use a single virtual element for each grain (polygons of any shape with n-nodes) with a significant reduction of the computational 70 burden, while the interconnecting layer is discretized with a fine mesh of triangular elements. Other advantages of VEM are: capability of using hanging nodes, that permits local refinements and coupling different degree elements; robustness to distortion of the elements; perfectly coupling with FEM elements; accuracy because the stiffness matrix is computed in precision machine; easy 75 implementation.
In this paper, low order virtual elements are used: the adoption of virtual elements of order one is particularly suitable for the homogenization procedure, as shown in [36,37,58]. Moreover, stress/strain is constant over the elements, but this approximation do not affect homogenization results. Furthermore, by 80 the adoption of the hanging nodes we subdivide the edges of the grains for generating a sufficiently fine mesh in the interphase part.
Exploiting FSHP combined with VEM, several parametric analyses have be performed to characterize overall mechanical parameters of polycrystalline composites with thin interphase layers and to investigate their sensitivity to the contrast 85 [59,60]. Attention is paid both at the identification of the RVE size and at the evaluation of the overall mechanical properties, that vary in relation with the contrast of the material.
The outline of the paper is as follows. In Section 2 the statistical procedure is recalled and specialized to the case at hand of polycrystals with thin inter-90 granular layers. Section 3 is devoted to the basic assumption of the first order virtual element formulation. Section 4 provides results of a set of parametric analyses on cermet-like polycrystals with intergranular phases stiffer of softer than grains. In order to check the capabilities of the proposed procedure, a successful comparison between the results obtained with VEM and those with a 95 very fine FEM mesh is performed. In Section 5 some final remarks are presented, highlighting the expected advantages of the proposed approach with respect to more standard approaches.

Fast Statistical Homogenization Procedure
The main features of the so-called Fast Statistical Homogenization Procedure 100 (FSHP), developed by the authors in [36,37], are detailed in this Section. First, the attention is focused on two types of composites, characterized either by circular inclusions embedded in a base matrix or by polycrystalline material with thin interfaces, for which a linear elastic constitutive behaviour is assumed at both the microscopic and the macroscopic level. Then, the key ideas of the 105 statistical homogenization procedure are briefly recalled and specialized for the materials at hand.

Computational homogenization
We take into account a two-dimensional continuum and describe, at the microscopic scale, the heterogeneous composite as a two-phase material. We 110 both investigate the case of composites made of circular inclusions of diameter d, randomly distributed in a matrix ( Fig. 1(a)), and the case of polycrystals (considered as inclusions) bounded by thin interfaces with average dimension d ( Fig. 1(b)). The former case is representative of ceramic/metal matrix composites or also concrete or rocks, while the latter one is representative of cermets.

115
In the present instance of non periodic composite materials and in view of the statistical homogenization procedure, it is useful to introduce a scale parameter δ = L/d * defined as the ratio between the edge of a square test window L, and the characteristic dimension d * . We furthermore define the material contrast 120 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 this case as low contrast materials, that are suitable to properly represent porous media [61]. On the other hand, when c > 1 inclusions are stiffer than the matrix and we refer to 125 high contrast materials [19].
In order to perform homogenization, we describe the material at two scales of interest: the microscopic and the macroscopic levels. At the microscopic level the  heterogeneous material is represented in detail, accounting for each constituent in terms of geometry and constitutive behaviour. At the macroscopic level the composite material is ideally replaced by an equivalent material whose global behaviour is representative of the actual heterogeneous material. The governing equations are formally the same as those defined at the microscopic level, except 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 micro-scale, while upper case letters to the macro-scale.
We refer to a linearised two-dimensional framework. At the lower level each material phase is characterized by linear elastic isotropic behaviour with the stress-strain relations written as: where ε ij and σ ij (i, j, k = 1, 2) are the components of the micro-strain and micro-stress tensors, λ and µ are the Lamé constants, and δ ij (i, j = 1, 2) is the Kronecker symbol. At the macroscopic level, the general anisotropic stressstrain relations, read: where E ij , Σ ij (i, j = 1, 2) are the components of the macro-strain and macrostress tensors and C ijhk (i, j, h, k = 1, 2) are the homogenized moduli, i.e the components of the macroscopic elastic tensor obtained via a homogenization procedure based on the Hill's macro-homogeneity condition [32]:  to ∂B δ (i, j = 1, 2).

Statistical Homogenization
FSHP is based on the statistical homogenization procedure previously developed in [19] and briefly described in Section 1. The proposed homogenization procedure is conceived for evaluating both the homogenized elastic parameters of a non-periodic heterogeneous material and for identify the Representative Volume Element (RVE), that in the absence of a repetitive micro-structure is 140 not known a priori.
According to the approach presented in [19], as well as in [21,11], the presented procedure requires the statistical definition of a number of realizations called Statistical Volume Elements (SVEs), representing the micro-structure, sampled in a Monte Carlo sense, which allows for determining series of scale-145 dependent upper and lower bounds for the overall elastic moduli and to approach the RVE size, δ RV E , using a statistical stopping criterion, which is based on the variation of the average elastic moduli.
All steps of the homogenization procedure are completely integrated in the FSHP and they are described below.

150
Step 1 Input: in the case of disk-shaped inclusion, we set set the nominal volume fraction of the medium (ρ ≤ 40%) as T ol, based on data dispersion, as defined in above. 165 Step 2 Input: initialize the window size, L = L 0 , and number of simulations, Step 3 Realizations: for the case of circular inclusion the procedure automatically determine for each window size the number of inclusions (Poisson random variable) via simulations exploit Knuth's algorithm (see [59,36]). PolyMesher developed by [62]. Then all edges are shifted inwards by s/2 175 in order to obtain the final geometry of all realizations B δ . In Fig. 2 a schematic of the procedure is depicted for a generic grain: the initial polygon has vertices P i ; additional pointsP i are identified as a result of the shifting procedure, so that the interface area (in light blue) and the grain area (in light red) are defined.  In all cases, each realization is supposed to be independent from any previ-   Step 4 Generate/Solve: for each SVE, generate the relative mesh and solve both the Dirichlet and Neumann (Eq. (4)) BVP, and compute the homogenized constitutive parameters.
Step 5 Compute: evaluate the average bulk modulus, K δ , the relative standard deviation σ( K δ ) and variation coefficient CV ( K δ ). Then com- 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, T ol. Repeat Steps 3-4 until N i < N lim .
Step 6 Checking: if the number of realizations necessary for ensuring the require-195 ment 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 whole procedure has been schematized in the flow-chart in Fig. 4. The statistical convergence criterion adopted is based on a 95% confidence level of the Normal Standard distribution, which provides the number N of realizations 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 210 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, T ol, are chosen in order to define a narrow con-215 fidence interval for the average and to obtain a reliable convergence criterion. The adopted statistical criterion allows us to detect the RVE size also when the Dirichlet and Neumann solutions do not tend to the same value. The values of the tolerance are assumed as a function of the data dispersion [19].
For all realizations of the micro-structure the Virtual Element Method (VEM) 220 has been adopted as a numeric tool, that permits to adopt single polygonal element for the inclusions without meshing with consequently high reduction of the computational burden with respect to finite elements. Moreover, the adoption of the VEM allows us to create a refined mesh in the interface without meshing the granular elements thanks to the use of so-called hanging nodes 225 (Fig. 5). In the zones of interfaces a triangular virtual element mesh have been adopted, determined using a code for generating random mesh of Delaunay type (Triangle, [63]). As demonstrated in [64], three nodes virtual elements behaves just like 3-noded triangular finite elements.
The computational strategies adopted are aimed at making the statistical ho-

Virtual Element Framework
In this section we briefly recall the governing equation of the linear elastostatic problem and the related weak formulation, mandatory starting point to describe 235 the construction of virtual element. We restrict our investigation to virtual elements of degree 1 [38,39] also referred as lower order virtual elements. The Voigt notation is adopted since is more suited to easily develop the VEM formulation.

The linear elastic problem
We consider a two-dimensional domain, Ω, with Γ being its boundary. The 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 the continuous bilinear form a(·, ·) : V × V → R, in witch R is the set of 250 the real numbers, reads: in witch C is the 3 × 3 elastic tensor, and the linear functional F(·) : V → R, reads:

Virtual element formulation
In this subsection we introduce the virtual element discretization used in the homogenization procedure. In order to approximate the solution of the problem (6) we consider a decomposition T h of the domain Ω into non overlapping 255 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. 6(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, 260 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 (6), as in the finite element methodology, in the following way: The discrete virtual element space, V h , is: where the local space, 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 h|E then is: We can observe that, in contrast to the standard finite element approach, the local space V h|E is not fully explicit. Moreover v h is a polynomial of degree 1 on each edge e of E and globally continuous on ∂E. The problem (6) 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. In this work the body force is null, so no details for its imple-mentation are reported. The discrete bilinear form is constructed element by element as By the above definition of the local space (11), the following important obser-265 vations can be made: -the functions v h ∈ V h|E are explicitly known on ∂E (linear functions); -the functions v h ∈ V h|E are not explicitly known inside the element E; The related degrees of freedom for the space V h|E are 2n e point-wise values

Projection operator and construction of the stiffness matrix
In accordance to [65] we define the projector operator of the strain as: More specifically, the projector operator respects the following orthogonality condition: As in the finite element method, we define the vectors v h ∈ V h|E andε ∈ P 0 (E) 2×2 sym using the degrees of freedom v h ∈ R 2ne , ε ∈ R 3 : The matrices N u and N V are where φ i is the standard i−th shape function. The projection operator in matrix is: with Π ∈ R 3×m .
By putting Eqs. (17), (19) into the Eq. (16) we obtain: Integrating by part the right hand side we obtain: where N E contains the components of the outward normal n = {n 1 , n 2 } T , i.e.
After some algebra: Finally the local projector operator Π is: where: We can observe that the matrix G ∈ R 3×3 is computed by knowing the area of the elements | E |: where | E | is computed using the Gauss-Green formula: where ∂E + is the boundary of the element oriented anticlockwise and the matrix where the boundary integrals in the matrix B are computed using the degrees of freedom of the element. In particular, the boundary integrals are computed resorting to the Gauss-Lobatto quadrature rule with k + 1 points, since this choice has turned to be adequate. In the case of VEM k = 1 we can obtain: where l i =| 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. Moreover if The main difference between Finite Element Method and Virtual Element Method regards the approximation of the bilinear form of Eq. (14), that for the VEM is as follow where the first term of the right hand side is the consistent term and the second term is the stabilization term. The stabilization part is not presented in Finite Element Method and it is one of the peculiarity of the VEM. For computing the consistent term M we use Eq. (19) and Eq. (7) where the stiffness matrix M is defined as Eq. (32) has been rearranged in compact form as follow adopting a constant elastic tensor Concerning the stabilization term, we adopt the choice introduced in [38,39].
However other authors have proposed a different stabilization terms [66]. Some preliminaries are required before introducing the stabilization term. In particular we introduce the space of the scaled monomials p ∈ (P k (E)) 2 since In the case of VEM of degree k = 1 the space of the scaled monomials is are the coordinates of the centre of E and h E its diameter. The stabilization term is defined as follow where τ ∈ R is a coefficient defined by the user that can be set equal to 1/2 for linear elasticity [65] and tr (·) denotes the trace operator. The matrix D ∈ R m×6 collects the values of the polynomials p i at the degrees of freedom on E (nodes of the element in the case of degree k = 1). An accurate inspection of Eq.

280
(35) reveals that the stability term is basically a rough approximation of the internal energy associated with the difference between VEM shape function and its projection, while τ tr(M ) scales the term with respect to the consistency part, which is the properly used for providing convergence of the method [38].

Local load vector
Regarding the approximation of right hand side of the equation (13), the same assumption proposed in [65] is used, i.e.
so that the body load over the element E is subdivided over all nodes υ i of the polygonal element.

Numerical Results
In this section, the results of parametric analyses performed using the described computational FSHP in conjunction with VEM are provided. We restrict our 290 investigation to polycrystalline material with thin interfaces, since the other type of composite has been already extensively studied in [67,58,68].

Comparison between VEM and FEM discretization
As a preliminary test, we consider an ideal material that mimics a cermet-like composite with grains and thin interfaces and investigate, for a given window The components of the homogenized elastic tensor C ijhk are evaluated using both Dirichlet (apex D) and Neumann (apex N) boundary conditions. The structure of the computed homogenized stiffness matrix is equal to: More specifically, the results obtained using with virtual elements are compared 310 to those obtained using finite elements, see

FSHP applied to materials A and B
In this example we take into account both materials A (stiff grains) and B (soft grains) with elastic properties reported in Tab. 3, where also the contrast c is 325 reported, and the average dimension of grains and the hard-core have been fixed equal to d = 40 µm and s = 2 µm, respectively.
As an example, in Fig. 8    The convergence trend for the different materials depends on the different dispersion of results, as shown in Fig. 10 and Fig. 11, where the Coefficient of Variation, CV (K), is also plotted for several values of window size L for Dirichlet (blue solid line) and Neumann (red solid line) BVPs.

340
In this paper we present an efficient extension of the Fast Statistical Homogenization (FSHP), previously developed by the authors for the homogenization of particle composites in [36,37,69], to account for the peculiar geometry of the polycrystal and composites with grains and thin interfaces. The Virtual Element Method is used as numerical tool to solve the partial differential equa-345 tions governing the elastic problem. Starting from a random polygonal mesh, a novel automatic procedure to insert thin interfaces between grains is presented.
Each grain is discretized using a single VEM, while the interfaces zones present a dense triangular mesh (exploiting Delaunay triangulation). Clear advantages in terms of computational burden are found with respect to the Finite Element for which large window domains are required, as reported in [30,19,34,9].
The results shown for a set of ideal materials, characterized by different mismatches between elastic moduli of grains and interfaces (material contrast), order to better describe the behaviour of micro-structured material, is one of 360 the future developments. Moreover, we will upgrade the procedure to take into account possible anisotropic behaviours, since this type of materials are gener-ally characterized by anisotropic properties, that varies grain-by-grain. Another interesting application and future development will be applying the procedure for the design of new materials with prescribed mechanical characteristics.