A coupled virtual element-interface model for analysis of fracture propagation in polycrystalline composites

Contents lists available at ScienceDirect Comput. Methods Appl. Mech. Engrg. journal homepage: www.elsevier.com/locate/cma A coupled virtual element-interface model for analysis of fracture propagation in polycrystalline composites Cristina Gatta, Marco Pingaro, Daniela Addessi, Patrizia Trovalusci∗ Department of Structural and Geotechnical Engineering, Sapienza University of Rome, Rome, Italy ARTICLE INFO Keywords: Random microstructure Fracture Polycrystalline composite Interface finite element Virtual element method ABSTRACT This paper proposes a coupled virtual element-interface finite element model for the analysis of the fracture propagation in polycrystalline composites with random microstructure. The key idea is to discretize each crystal, also referred to as grain, with a single low order virtual element with elastic constitutive response, and describe the interaction between grains by means of damaging and frictional zero-thickness interface finite elements. Thus, the typical intergranular crack growth is modeled by avoiding refined finite element grain discretizations with relevant computational cost saving. Results of numerical simulations are presented and discussed. First, some benchmarks show the reliability of the proposed modeling strategy. Then, the response of Alumina/Zirconia representative volume elements, whose size is selected on the basis of results of a statistical homogenization procedure tailored for random composites, is investigated by analyzing the effect of the variation of the metallic phase volume fraction and the shape of grains composing the microstructure. 1. Introduction Ceramic matrix composites (CMCs) are a special class of materials, comprising ceramic particles or short/long fibers dispersed within a ceramic matrix, which can also include polycrystalline or layered structures. These composites are engineered to overcome the limitations commonly observed in conventional ceramics used in various technical applications, particularly concerning their brittle-like fracture under mechanical or thermal stresses [1]. The inclusion of ceramic particles or fibers enhances fracture toughness, potentially transitioning the brittle response to a more ductile fracture behavior and improving the material resistance to thermal shock. Moreover, CMCs retain the advantageous properties of the ceramic matrix, such as high strength and Young’s modulus. In addition to their increased fracture toughness and resistance to thermal shock, CMCs also exhibit excellent properties of corrosion resistance and high temperature endurance. Their capability to maintain mechanical properties at extreme temperatures makes them ideal for high-temperature and high-pressure environments. Furthermore, CMCs can be designed to have a lighter specific weight compared to metallic materials, offering significant advantages in the aerospace and automotive applications, where weight is a critical factor [2–5]. The corrosion resistance is important, especially in aggressive environments such as marine or chemically reactive ones, where long-term durability is crucial. Additionally, the dimensional stability of CMCs, combined with their wear and fatigue resistance, makes them suitable for high dynamic load applications, such as rotating components in high-performance machinery. Hence, the applications of CMCs are manifold and encompass critical areas such as heat shield ∗Corresponding author. E-mail addresses: [email protected] (C. Gatta), [email protected] (M. Pingaro), [email protected] (D. Addessi), [email protected] (P. Trovalusci). https://doi.org/10.1016/j.cma.2024.117383 Computer Methods in Applied Mechanics and Engineering 432 (2024) 117383 Available online 28 September 2024 0045-7825/Published by Elsevier B.V. This is an open access article under the CC BY license ( http://creativecommons.org/licenses/by/4.0/ ). C. Gatta et al. systems for space vehicles, brake disks, slide bearings, components for high-temperature gas turbines, cutting tools, and burner components [6–9]. An outstanding prototype of polycrystalline CMC is undeniably the Alumina and phase-stabilized Zirconia composite Al2O3/ZrO2 (with varying volume contents), as discussed in [10]. This composite merges the favorable features of Alumina, such as high hardness and low susceptibility to aging, and Zirconia, namely high fracture toughness and resistance to subcritical crack growth [5,11]. The Al2O3/ZrO2composite is characterized by a complex internal arrangement composed of Alumina and Zirconia grains with dimensions varying from 0.2 μm to 2 μm. A cluster of randomly dispersed grains forms a polycrystalline structure. Its mechanical characterization has garnered considerable attention within the scientific community, as evidenced by [12–18], with the primary objective of designing optimized materials capable of meeting high-tech demands. Relying on the above considerations, this paper focuses on the nonlinear response of Alumina/Zirconia Representative Volume Elements (RVEs), considering the influence of different volume content of the two components and accounting for the effect of the microstructure randomness. Specifically, we propose an efficient and reliable nonlinear model that merges Virtual Elements (VEs) and Interface Finite Elements (IFEs), aiming to reduce computational burden although preserving accuracy. The Virtual Element Method (VEM) is a numerical procedure [19,20] developed in recent years as a viable alternative to the Finite Element Method (FEM), for a broad spectrum of mechanical applications [21–26]. The method has recently been enhanced [27,28] and extended to perform large displacement analysis within the corotational framework [29,30]. Its peculiar features, in terms of high flexibility in the number of nodes and shape of elements, make it suitable for modeling materials with polycrystalline microstructures [31,32], as each grain can be discretized with just one virtual element with generic shape. However, it should be mentioned that the literature also proposed other methods capable of managing polygonal elements to mimic the natural shape of each grain, such as Polygonal Finite Elements [33–37], Voronoi cell Finite Elements [38–41], and Trefftz–Lekhnitskii Grains [42]. Interface elements have been employed in many engineering fields to study a variety of problems, including adhesive bonding, mechanical contacts and fracture propagation. Overall, interface mechanics establishes the link between materials (or structures) bonded together by cohesive forces with various nature. In the framework of the FEM, zero-thickness interface elements have been developed, based on one-, twoand three-dimensional formulations, considering both material and geometric nonlinearities [43–45]. A comprehensive review of interesting interface models that consider the combined effect of damage, unilateral contact and friction can be found in [46]. Recently, virtual elements and interface elements have been coupled to analyze the processes of nucleation and evolution of cohesive fracture in two-dimensional bodies [47] and, also, to simulate cracking of cement-based composites [48]. In this work, we employ damage-friction IFEs to reproduce the fracture propagation phenomenon typically occurring at the intergranular zones in CMC materials beyond the elastic regime. This study employs a two-dimensional (2D) formulation. Due to their geometric characteristics and boundary conditions, a number of significant cases can be effectively captured using 2D formulations, relying on well-considered assumptions. While threedimensional (3D) models offer a more precise and detailed depiction of material behavior, they are often computationally demanding and complex to develop. In contrast, adopting a 2D modeling approach not only simplifies the interpretation of numerical results but also enables the execution of a large number of analyses, providing deeper insights into the material’s nonlinear response [48–51]. Overall, 2D models are a good choice for qualitative intuitions, but their use should be done with caution, ensuring that appropriate geometric features and boundary conditions are valid. It is widely recognized that identification of the RVE size is well established in periodicity-based homogenization techniques resorting to single cell concept, while it is still an opened challenge for random media [52]. These procedures are employed within limit processes involving several finite-scale continuous descriptions (Statistical Volume Elements), relative to the microstructural length scale. The solution of series of Dirichlet and Neumann boundary value problems (BVPs) at several mesoscales deduced from Hill–Mandel type macrohomogeneity condition, also valid for non-periodic and non classical media [53], provides two hierarchies of bounds for the material properties, as well as the microstructural minimal size of the RVE for performing homogenization [52,54,55]. To address this issue, we employ the Fast Statistical Homogenization Procedure (FSHP) in conjunction with the Virtual Element Method, developed by some of the authors [56,57]. This allows to rapidly determine the RVE for the CMC material operating within the framework of a first-order computational homogenization scheme. Both the equivalent elastic moduli and the characteristic size of the RVE are determined by establishing bounds on the effective response. These bounds are obtained by solving boundary value problems with either Dirichlet or Neumann boundary conditions, as discussed in [58,59], following approaches similar to those presented in [60–62] within the context of micropolar continua. The paper is organized as follows. Section 2recalls the basics of the statistical homogenization applied to the CMC material to determine the RVE size to be used in the following nonlinear analyses. Section 3is dedicated to the model description, giving details on the virtual element and interface formulations adopted. Section 4illustrates some benchmark tests finalized to validate the proposed modeling strategy. Section 5investigates the response of Alumina/Zirconia RVEs focusing on the effect of the variation of microstructure and volume fraction of the two materials. Finally, concluding remarks are given in Section 6. 2. From statistical homogenization to fast statistical homogenization procedure The main aim of this section is to recall the basics of the so-called Fast Statistical Homogenization Procedure (FSHP), already developed by some of the authors in [56,57] for particular topology, i.e. composites made of circular inclusions, representative of fibers, randomly dispersed in a second phase, called matrix. The peculiar structure of the CMC composites, characterized by random geometry of the particles besides random positions, required an ad-hoc procedure for generation of the Statistical Volume Element Computer Methods in Applied Mechanics and Engineering 432 (2024) 117383 2 C. Gatta et al. Fig. 1. Neumann boundary conditions applied on the window 𝛿. Source: Modified by [32]. (SVE), as developed in [32]. Here, FSHP inspired by the pioneering works of [54,61] is adopted to determine the RVE for CMC material. In the following, taking into account the symmetries of the tensors, in 2D framework we represent the second order tensors as 3-component vectors and the fourth order tensors as 3x3 matrices. 2.1. First order homogenization A classical energy-based computational homogenization approach [63,64] has been adopted as a step of the Statistical Homogenization Procedure with the aim of estimating the components of the overall elastic matrix associated to random realizations of the CMC Alumina/Zirconia. To perform the homogenization, we describe the material at two scales of interest: the microscopic and macroscopic levels. At the microscopic level, the heterogeneous material is represented in detail, accounting for all the geometric and constitutive properties of the constituents. At the macroscopic level, the composite material is ideally replaced by an equivalent medium, whose global response is representative of that of the actual heterogeneous material. All the governing equations are known at this level and 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 the result of the homogenization procedure. In the following, lower case letters are always related to the micro-scale, while upper case letters to the macro-scale. In view of the statistical homogenization procedure, it is useful to introduce a scale parameter 𝛿=𝐿∕𝑑defined, at the microscopic scale, as the ratio between the edge length, 𝐿, of a square test window, and the characteristic dimension, 𝑑, of a grain. Referring to a 2D framework, at the lower level each material phase is characterized by linear elastic isotropic behavior with the stress–strain relations written as: 𝝈=C𝜺(1) where 𝜺and 𝝈are micro-strain and micro-stress vectors, and Cis the elastic constitutive matrix at the microscale. At the macroscopic level, the general anisotropic stress–strain relation reads: 𝜮=C𝐄(2) where 𝐄,𝜮are the macro-strain and macro-stress vectors and Cis the homogenized material stiffness matrix that contains the homogenized moduli : C=⎡⎢⎢⎢⎣ C1111 C1122 C1112 C2211 C2222 C2212 C1211 C1222 C1212⎤⎥⎥⎥⎦ (3) i.e. the components of the macroscopic elastic matrix obtained via a homogenization procedure based on the Hill–Mandel macro-homogeneity condition [65]: 𝐄𝑇𝜮=1 𝐴𝛿∫𝛿 𝜺𝑇𝝈𝑑𝐴 (4) which establishes an equivalence of the average internal work over 𝛿(test window), occupying a region of area 𝐴𝛿=∫𝛿𝑑𝐴, and the mechanical internal work density of the macroscale model, expressed in terms of homogenized stress and strain measures. The homogenization procedure is based on the solution of properly defined boundary value problems at the microscopic level with Neumann (Fig. 1) and Dirichlet (Fig. 2) boundary conditions applied on the boundary 𝜕𝛿, directly deriving from the fulfillment of the macro-homogeneity condition in Eq. (4). Computer Methods in Applied Mechanics and Engineering 432 (2024) 117383 3 C. Gatta et al. Fig. 2. Dirichlet boundary conditions applied on the window 𝛿. Source: Modified by [32]. 2.2. Fast statistical homogenization The FSHP is based on the statistical homogenization procedure proposed in [66], then developed for micropolar continua in [61] and recently automatized in [56,57]. The procedure is conceived both for evaluating the homogenized elastic parameters of a non-periodic heterogeneous material, and identifying the RVE, which is not known a priori in the absence of a repetitive microstructure. According to the approach presented in [52,61,67], the hypotheses of statistical homogeneity and isotropy, combined with the mean-ergodicity of the microstructure, are assumed. In this framework, the presented procedure requires the statistical definition of a number of realizations called Statistical Volume Elements (SVEs), representing the microstructure, sampled in a Monte Carlo sense, which allows for determining series of scale-dependent upper and lower bounds for the overall elastic moduli and to approach the RVE size, corresponding to 𝛿𝑅𝑉 𝐸 , using a statistical stopping criterion based on the variation of the average elastic moduli. All steps of the homogenization procedure are completely integrated in the FSHP. They are described below and schematically depicted in the flow-chart in Fig. 3. Step 1 Input: set the average size of the grains 𝑑and define the dimensionless scale factor 𝛿=𝐿∕𝑑. Fix the mechanical parameters of each phase, i.e. Young’s modulus and Poisson coefficient of each phase, 𝐸𝑖and 𝜈𝑖(𝑖=1,2). Set the minimum number of simulations for convergence, 𝑁𝑙𝑖𝑚, and a tolerance parameter, 𝑇 𝑜𝑙, based on data dispersion. Step 2 Input: initialize the window size, 𝐿=𝐿0, and number of simulations, 𝑁=𝑁0. Step 3 Realizations: generate a random polygonal mesh with average size of grains 𝑑using the MATLAB®program PolyMesher developed by [68]. Each realization is supposed to be independent from any previous one. Based on volume fractions, mechanical parameters are randomly assigned to grains. In order to avoid abnormal boundary layers related to the artefact of generating the Voronoi tessellations from realizations of homogeneous random point fields created only within the considered test windows, we generate the realizations used in the numerical simulations by cutting out smaller windows. In Fig. 4, a realization obtained by FSHP for different windows sizes has been plotted, highlighting in green the cutting windows. Step 4 Generate/Solve: for each SVE, generate the relative mesh and solve both the Dirichlet and Neumann BVPs, and compute the homogenized constitutive parameters. Step 5 Compute: the homogenized bulk modulus K=((C1122 +C2211)∕2 + C1212)∕6 and evaluate the average bulk modulus, ⟨K⟩𝛿, the relative standard deviation (⟨K⟩𝛿)and variation coefficient 𝐶𝑉 (⟨K⟩𝛿) defined as the ratio of the standard deviation to the mean 𝜇 (𝐶𝑉 =∕𝜇). Then, compute 𝑁𝑖=(1.96 (⟨K⟩𝛿) ⟨K⟩𝛿𝑇 𝑜𝑙 )2 (5) 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, 𝑇 𝑜𝑙. Repeat Steps 3–4 until 𝑁𝑖⩽𝑁𝑙𝑖𝑚. We focus our attention on the bulk modulus as a representative material parameter because the homogenized elastic tensor has mainly found to be isotropic [31,56]. Step 6 Checking: if the number of realizations necessary to ensure the requirement at Step 5 is small enough, stop the procedure. We choose as the needed number of realizations the largest unfavorable number between those obtained by solving BVPs of Neumann or Dirichlet. Otherwise, choose an increased value of 𝛿and go to Step 3. Computer Methods in Applied Mechanics and Engineering 432 (2024) 117383 4 C. Gatta et al. Fig. 3. Flow-chart of the fast statistical homogenization procedure. Source: Adapted from [69]. Fig. 4. Example of realizations obtained by FSHP for different window size 𝛿=𝐿∕𝑑with in green highlighted the cutting windows. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) The fulfillment 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 set. The effective homogenized elastic moduli can be determined as the arithmetic mean value between the Dirichlet (upper) and Neumann (lower) bounds at the convergence window (RVE). The statistical convergence criterion adopted is based on a 95% confidence level of the Normal Standard distribution, which provides the number 𝑁of realizations at which it 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 Computer Methods in Applied Mechanics and Engineering 432 (2024) 117383 5 C. Gatta et al. Fig. 5. Centroidal Voronoi tessellation microstructure generated via PolyMesher [68] (left) and Virtual Element extract to the mesh (right). Fig. 6. Examples of realizations with different content of Alumina and Zirconia for window dimension 𝛿= 9 (Alumina and Zirconia grains are depicted in green and blue color, respectively). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) to reaching the minimum window size 𝛿𝑅𝑉 𝐸 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, 𝑁𝑙𝑖𝑚, and the tolerance parameter, 𝑇 𝑜𝑙, are chosen in order to define a narrow confidence interval for the average and 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 [61]. As already shown in [31], the geometrical features of a polycrystalline material are particularly suitable for using the Virtual Element Method as a numerical tool for solving BVPs at the microscopic scale. In this case, indeed, the microstructure can be satisfactorily modeled with a randomly generated centroidal Voronoi tessellation and it can be directly used as a computational mesh. As is well known, in fact, the recent numerical tool of the VEM permits to use single polygonal element for the grains (Fig. 5), avoiding internal meshing, with consequent significant reduction of the computational burden with respect to finite elements. Note that at this stage, for the sake of simplicity, the average size 𝑑of grains is kept constant, but a straightforward modification to account for different sizes is possible. The computational strategies adopted are aimed at making the statistical homogenization process as efficient as possible for solving a series (hundreds) of boundary value problems, required by the statistical homogenization procedure and to rapidly converge to the RVE solution. The capability of the VEM of delivering reliable estimations of overall elastic moduli, despite the use of very coarse meshes, has been already assessed in [31,56]. 2.3. Identification of Alumina/Zirconia RVE The FSHP described is applied to identify the Alumina/Zirconia CMC material RVE, taking into account the influence of randomness of the microstructural topology. Linear elastic isotropic material is assumed for each phase, setting Young’s moduli and Poisson ratios of Zirconia (Z) and Aluminia (A) equal to 𝐸Z= 138 GPa,𝐸A= 345 GPa,𝜈Z= 0.28 and 𝜈A= 0.25, respectively [17]. Four different Al2O3/ZrO2composites are considered. These are all characterized by the average grain size 𝑑= 0.4 μm, but different volume fraction of Alumina, 𝑝, which varies in the range [20%÷80%].Fig. 6 shows an example of realizations corresponding to the four types of material examined, when the dimensionless scale factor is 𝛿= 9. The results of the homogenization procedure prove that the overall behavior of the composite material does not significantly differ from that of an equivalent isotropic material. Moreover, the convergence trend for the different materials depends on the different dispersion of results, as shown in Figs. 7(a)–8(a)–9(a)–10(a), where the Coefficient of Variation, CV, is plotted versus 𝛿for Dirichlet (blue solid line) and Neumann (green solid line) solutions. Computer Methods in Applied Mechanics and Engineering 432 (2024) 117383 6 C. Gatta et al. Fig. 7. Al2O3𝑝= 20%. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) Fig. 8. Al2O3𝑝= 40%. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) Applying the convergence criterion in Eq. (5), it emerges that results reach the convergence values at 𝛿= 10, which defines the dimension of the RVE. This is confirmed by the number of simulations needed to achieve convergence by varying the window size 𝛿(Figs. 7(b)–8(b)–9(b)–10(b)). This result is adopted in the following nonlinear analyses. 3. Formulation of the coupled virtual element-interface element model In this section we describe the nonlinear coupled virtual element-interface finite element model proposed. First, the governing equations of the linear elastostatic problem and its approximation using the Virtual Element Method are recalled, then, the attention is dedicated to the interface formulation. The model is developed in the 2D framework under plane-stress assumptions and the hypothesis of small displacements and strains. Computer Methods in Applied Mechanics and Engineering 432 (2024) 117383 7 C. Gatta et al. Fig. 9. Al2O3𝑝= 60%. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) Fig. 10. Al2O3𝑝= 80%. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) 3.1. Virtual element formulation We consider a 2D domain 𝛺 ⊆ R2with piecewise smooth boundary 𝜕𝛺, as depicted in Fig. 11. The body is subjected to the volume force, represented by the vector 𝐟∈(𝐿2(𝛺))2,𝐟= {𝑓1𝑓2}𝑇(where 𝐿2(𝛺)is the standard Lebesgue space), and given boundary conditions, under the hypothesis of small strains. For the sake of simplicity, we consider homogeneous Dirichlet boundary conditions. We introduce the space of admissible displacement field, 𝐕∶= (𝐻1 0(𝛺))2(with 𝐻being the Sobolev space), represented by the vector 𝛿𝐮. Furthermore, we represent the strain as the vector 𝜺(𝐮)associated to the displacement field 𝐮= {𝑢1𝑢2}𝑇as: 𝜺(𝐮) = 𝐒 𝐮 with 𝐒=⎡⎢⎢⎣ 𝜕𝑥10 0𝜕𝑥2 𝜕𝑥2𝜕𝑥1⎤⎥⎥⎦ (6) where 𝜕(⋅)denotes the partial derivative operator with respect to the (⋅)-coordinate. Computer Methods in Applied Mechanics and Engineering 432 (2024) 117383 8 C. Gatta et al. Fig. 11. Two-dimensional domain (left) and its tessellation in non overlapping polygons (right). By using the Principle of Virtual Work it is possible to define the weak form of the linear elasto-static problem, which reads: {Find 𝐮∈𝐕such that : 𝑎(𝛿𝐮,𝐮) = (𝛿𝐮) ∀ 𝛿𝐮∈𝐕(7) where 𝑎(𝛿𝐮,𝐮) = ∫𝛺 𝜺(𝛿𝐮)𝑇C𝜺(𝐮)𝑑𝛺 (8) is the continuous bi-linear form 𝑎(⋅,⋅) ∶ 𝐕×𝐕→R. Furthermore, the linear functional (⋅) ∶ 𝐕→R, reads: (𝛿𝐮) = ∫𝛺 𝛿𝐮𝑇𝐟𝑑𝛺 (9) In order to approximate the solution of the problem in Eq. (7), we consider a decomposition ℎof the domain 𝛺into non overlapping polygonal elements 𝐸(Fig. 11). In the following, we denote by 𝑒the straight edges of the mesh ℎfor all 𝑒∈𝜕𝐸. The symbol 𝑛𝑒represents the number of the vertexes of the polygon 𝐸. Let 𝑘be an integer ≥1and 𝑘(𝛺)be the space of polynomials, living on the set 𝛺 ⊆ R2, of degree less than or equal to 𝑘. The discrete virtual element space 𝐕ℎ⊆𝐕is directly defined element by element and results: 𝐕ℎ∶= {𝛿𝐮∈𝐕∶𝛿𝐮|𝐸∈𝐕ℎ|𝐸∀ℎ}(10) being 𝐕ℎ|𝐸∶= [𝑉ℎ|𝐸]2and the local space defined as follows: 𝑉ℎ|𝐸∶= {𝛿𝑢ℎ∈𝐻1(𝐸) ∩ 𝐶0(𝐸) ∶ 𝛥𝛿𝑢ℎ∈𝑘−2(𝐸), 𝛿𝑢ℎ|𝑒∈𝑘(𝑒) ∀ 𝑒∈𝜕𝐸}(11) 𝑘−2(𝐸)is {0} when 𝑘= 1. The space dimension 𝐕ℎ|𝐸is 𝑚= 2𝑘𝑛𝑒+𝑘(𝑘− 1). Particularly, 2𝑛𝑒degrees of freedom are located at the vertexes of the polygon 𝐸,2𝑛𝑒(𝑘− 1) degrees of freedom are located at the edges 𝑒(2(𝑘− 1) in each edge) and 𝑘(𝑘− 1) degrees of freedom are the moments defined as 1 ∣𝐸∣∫𝐸𝒒𝑇𝒖ℎ𝑑𝐸 for 𝒒∈(𝑘−2(𝐸))2. The definition of the local virtual element space in Eq. (11) is the classical one proposed in [20] and here adopted, but there are other possible choices [27,29]. Focusing on the definition of the local space 𝑉ℎ|𝐸, it is remarked that the functions 𝛿𝐮ℎare not explicitly known inside the element but only on the element boundary. Considering the space discretization in Eq. (10) the discrete form of the weak formulation of the problem in Eq. (7) results as follows: {Find 𝐮ℎ∈𝐕ℎsuch that: 𝑎ℎ(𝛿𝐮ℎ,𝐮ℎ) = (𝛿𝐮ℎ) ∀ 𝛿𝐮ℎ∈𝐕ℎ (12) where 𝑎ℎ(⋅,⋅) ∶ 𝐕ℎ×𝐕ℎ→Ris the discrete bi-linear form approximating the continuous form 𝑎(⋅,⋅). This is built element by element, resulting as: 𝑎ℎ(𝛿𝐮ℎ,𝐮ℎ) = ∑ 𝐸∈ℎ 𝑎𝐸 ℎ(𝛿𝐮ℎ,𝐮ℎ) ∀ 𝛿𝐮ℎ,𝐮ℎ∈𝐕ℎ(13) The term (𝛿𝐮ℎ)in Eq. (12) approximates the virtual work of the external loads. Due to the peculiar definition of the local space 𝐕ℎ, the evaluation of the stiffness matrix goes through the definition of the so-called projector operator 𝛱, which is: 𝛱∶𝐕ℎ|𝐸→𝑘−1(𝐸)2×2 𝑠𝑦𝑚 𝛿𝐮ℎ→𝛱(𝛿𝐮ℎ)(14) Computer Methods in Applied Mechanics and Engineering 432 (2024) 117383 9 C. Gatta et al. Fig. 21. Alumina/Zirconia RVEs with different percentages 𝑝of Alumina. Fig. 22. Tensile response of Alumina/Zirconia RVEs with different percentages of Alumina: load–displacement global curves. 5.2. Effect of microstructure The effect of the grain shape and distribution, strictly connected to the randomness characteristic of the microstructure, is here investigated considering the three specimens depicted in Fig. 25. These are characterized by the microstructures denoted as ‘G1’, ‘G2’ and ‘G3’, but same Alumina volume fraction 𝑝= 60%. Note that microstructure G1 corresponds to that studied in Section 5.1. Assuming the loading and boundary conditions described in detail in Section 4,Figs. 26 and 27 show the uni-axial tensile response of the specimens in terms of load–displacement global curves and deformed configurations, respectively. The curves exhibit similar peak loads and slope of the softening branches but, also, some discrepancies due to the occurring stress distributions that induce different localizations of the main crack. Finally, the investigation is moved to another loading scenario. The assumed loading and boundary conditions correspond to those detailed in Section 4for the shear test. These completely restrain the base of the specimens and impose an increasing horizontal displacement at all nodes of their top side. Computer Methods in Applied Mechanics and Engineering 432 (2024) 117383 16 C. Gatta et al. Fig. 23. Tensile response of Alumina/Zirconia RVEs with different percentages of Alumina: cracks at the end of the analyses (scale factor of the deformed configuration equal to 100). Fig. 24. Tensile response of Alumina and Zirconia RVEs: distributions of the damage variable, 𝐷, at the interfaces in correspondence of points A, B, C, D, E, F and G indicated in Fig. 22. The resulting global curves (Fig. 28) show very similar features with no sudden drops of the force. In fact, for all the tested microstructures, the fracture process starts at the bottom left corner, where the maximum tensile stresses at the interfaces are concentrated, and, then, gradually spreads around crossing the entire RVE and inducing its collapse (Fig. 29). 6. Conclusions This work aimed to propose a computationally efficient modeling approach for the analysis of the fracture propagation in polycrystalline composites with thin interfaces or grain boundaries. The model exploits in a combined way the advantages of the modern and appealing Virtual Element Method and the consolidate Finite Element Method. The key idea is to discretize each crystal with a single virtual element and model the interaction behavior between the particles by means of a mesh composed of interface finite elements, where the damaging and frictional phenomena are assumed to be concentrated. This allows to capture the typical failure mechanism due to the intergranular fractures avoiding the dense FE discretizations usually employed. Computer Methods in Applied Mechanics and Engineering 432 (2024) 117383 17 C. Gatta et al. Fig. 25. Alumina/Zirconia RVEs with percentage of Alumina equal to 𝑝= 60% and different microstructures G1, G2 and G3. Fig. 26. Tensile response of Alumina/Zirconia RVEs with percentage of Alumina equal to 𝑝= 60% and different microstructures: load–displacement global curves. Fig. 27. Tensile response of Alumina/Zirconia RVEs with percentage of alumina equal to 𝑝= 60% and different microstructures: cracks at the end of the analyses (scale factor of the deformed configuration equal to 100). The modeling strategy is actually applicable to many categories of composites, with random or periodic microstructure. In this work, the attention was dedicated to the study of the micromechanical response of CMC representative volume elements, with special focus on the paradigmatic example of Alumina/Zirconia composites. The main contributions and novelty of the paper can be listed as follows: •the coupled use of virtual elements and interface elements has been applied for the first time to the modeling of CMC materials in the nonlinear field to track intergranular fracture propagation; •the developed formulation allows for significant mesh flexibility with low computational cost, as the admissible virtual elements are polygons of any shape, possibly characterized by hanging nodes, that well mimic the shape of particles of polycrystalline composites; Computer Methods in Applied Mechanics and Engineering 432 (2024) 117383 18 C. Gatta et al. Fig. 28. Shear test on Alumina/Zirconia RVEs with percentage of Alumina equal to 𝑝= 60% and different microstructures: load–displacement global curves. Fig. 29. Shear test on Alumina/Zirconia RVEs with percentage of Alumina equal to 𝑝= 60% and different microstructures: cracks at the end of the analyses (scale factor of the deformed configuration equal to 10). •the Fast Statistical Homogenization Technique proposed in [32,56] has been applied to select in reasonable manner the RVE of the random media analyzed; •the campaign of comparative analyses performed, involving tensile and shear tests, proved the reliability of the proposed approach and allowed to trace comparisons also between VEM and FEM computations. Results derived from coarse meshes, where the interface between grains is modeled with a single interface FE and each crystal is represented with a single VE, well represent those obtained using denser discretizations for grains and interfaces, with a maximum computational saving of approximately 87% in terms of analysis time; •the numerical simulations performed on Alumina/Zirconia RVEs highlighted the effects of the main factors affecting the response of these composites, i.e. the variation of the volume fraction of the material components and that of the microstructure. It emerged that, as the Zirconia content increases, a more brittle tensile response occurs with the load–displacement global curves exhibiting more severe softening branches. As for the influence of the microstructure, the response curves exhibit similar peak loads and slope of the softening branches, though different cracking paths occur under tensile loads in the analyzed cases. To summarize, the model appears as an efficient numerical tool for analyzing response of composite materials characterized by different representative size of the microstructure. For instance, this could be embedded in a multiscale procedure for the analysis of polycrystalline materials with nanoscale size grains randomly distributed. In this framework, the constitutive response of the homogeneous material considered at the macroscopic level could be derived by the step-by-step detailed analysis of the RVE identified through the FSHP [32,56] and modeled with the VE-FE formulation presented. Moreover, the FSHP could be extended to the nonlinear range, so that to evaluate the mechanical properties of the structural-scale material with an a-priori RVE-based homogenization. The model can be further improved by adopting advanced virtual element formulations accounting for enhanced description of the strain field leading to a self-stabilized elements [27,29]. Moreover, higher order continua, such as Couple-Stress and Cosserat, could be adopted for the grain description, as these have proved to better represent the behavior of microstructured material with random particles having different shape and orientation [80]. Computer Methods in Applied Mechanics and Engineering 432 (2024) 117383 19 C. Gatta et al. Finally, the decrease of the number of degrees of freedom offered by the use of coarse VE-meshes makes the adopted modeling strategy attractive to model real case microstructures. CRediT authorship contribution statement Cristina Gatta: Writing – original draft, Software, Methodology, Investigation, Conceptualization. Marco Pingaro: Writing – original draft, Validation, Software, Methodology, Investigation, Conceptualization. Daniela Addessi: Supervision, Funding acquisition, Conceptualization. Patrizia Trovalusci: Supervision, Funding acquisition, Conceptualization. Declaration of competing interest The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: Patrizia Trovalusci reports financial support was provided by University of Rome La Sapienza. If there are other authors, they declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Data availability No data was used for the research described in the article. Acknowledgments The authors acknowledge the grant PNRR PE5-CHANGES-Spoke 7 (CUP: B53C22003780006) and the grant PNRR CN1-Spoke 6 (CUP: B83C22002940006). References [1] A. Okada, Automotive and industrial applications of structural ceramics in Japan, J. Eur. Ceram. Soc. 28 (5) (2008) 1097–1104. [2] X. Tian, C. Birk, C. Du, E.T. Ooi, Automatic micro-scale modelling and evaluation of effective properties of highly porous ceramic matrix materials using the scaled boundary finite element method, Comput. Methods Appl. Mech. Engrg. 419 (2024) 116596, http://dx.doi.org/10.1016/j.cma.2023.116596. [3] N.P. Padture, Advanced structural ceramics in aerospace propulsion, Nature Mater. 15 (8) (2016) 804–809, http://dx.doi.org/10.1038/nmat4687. [4] W.G. Fahrenholtz, E.J. Wuchina, W.E. Lee, Y. Zhou, Ultra-high temperature ceramics: Materials for extreme environment applications, vol. 9781118700785, 2014, pp. 1–441, http://dx.doi.org/10.1002/9781118700853, [5] T. Sadowski, L. Marsavina, Multiscale modelling of two-phase ceramic matrix composites, Comput. Mater. Sci. 50 (4) (2011) 1336–1346. [6] F. Christin, CMC materials for space and aeronautical applications, Ceram. Matrix Compos.: Fiber Reinf. Ceram. Appl. (2008) 327–351. [7] W. Krenkel, Ceramic Matrix Composites: Fiber Reinforced Ceramics and Their Applications, John Wiley & Sons, 2008. [8] F. Raether, Ceramic matrix compositesAn alternative for challenging construction tasks, Ceram. Appl. 1 (1) (2013) 45–49. [9] P. Spriet, CMC applications to gas turbines, Ceram. Matrix Compos.: Mater. Model. Technol. (2014) 591–608. [10] V. Naglieri, P. Palmero, L. Montanaro, J. Chevalier, Elaboration of alumina-zirconia composites: Role of the zirconia content on the microstructure and mechanical properties, Materials 6 (5) (2013) 2090–2102. [11] M. Boniecki, T. Sadowski, P. Gołębiewski, H. Węglarz, A. Piątkowska, M. Romaniec, K. Krzyżak, K. Łosiewicz, Mechanical properties of alumina/zirconia composites, Ceram. Int. 46 (1) (2020) 1033–1039, http://dx.doi.org/10.1016/j.ceramint.2019.09.068. [12] O.Y. Zadorozhnaya, T. Khabas, O. Tiunova, S. Malykhin, Effect of grain size and amount of zirconia on the physical and mechanical properties and the wear resistance of zirconia-toughened alumina, Ceram. Int. 46 (7) (2020) 9263–9270, http://dx.doi.org/10.1016/j.ceramint.2019.12.180. [13] T. Waqar, S.S. Akhtar, A.F.M. Arif, A.S. Hakeem, Design and development of ceramic-based composites with tailored properties for cutting tool inserts, Ceram. Int. 44 (18) (2018) 22421–22431, http://dx.doi.org/10.1016/j.ceramint.2018.09.009. [14] T. Norfauzi, A. Hadzley, U. Azlan, A. Afuza, M. Faiz, M. Naim, Fabrication and machining performance of ceramic cutting tool based on the Al2O3-ZrO2-Cr2O3compositions, J. Mater. Res. Technol. 8 (6) (2019) 5114–5123, http://dx.doi.org/10.1016/j.jmrt.2019.08.034. [15] A. Afuza, A. Hadzley, T. Norfauzi, A. Umar, M. Faiz, F. Naim, K. Thongkaew, Analysis of particles size distribution on the agglomeration and shrinkage of alumina-zirconia compacts, Int. J. Nanoelectron. Mater. 13 (Special Issue ISSTE2019) (2020) 277–284. [16] W. Zhao, J. Chang, Q. Wei, J. Wu, C. Ye, Effect of sodium silicate solution combined yttrium oxide stabilized zirconia nanopowders on the properties of alumina ceramics fabricated by binder jetting additive manufacturing, J. Mater. Process. Technol. 330 (2024) http://dx.doi.org/10.1016/j.jmatprotec. 2024.118454. [17] T. Sadowski, B. Pankowski, Numerical modelling of two-phase ceramic composite response under uniaxial loading, Compos. Struct. 143 (2016) 388–394. [18] T. Sadowski, K. Łosiewicz, M. Boniecki, M. Szutkowska, Assessment of mechanical properties by nano-and microindentation of alumina/zirconia composites, Mater. Today: Proc. 45 (2021) 4196–4201. [19] L. Beirão Da Veiga, F. Brezzi, A. Cangiani, G. Manzini, L. Marini, A. Russo, Basic principles of virtual element methods, Math. Models Methods Appl. Sci. 23 (1) (2013) 199–214, http://dx.doi.org/10.1142/S0218202512500492. [20] L. Beirão Da Veiga, F. Brezzi, L. Marini, Virtual elements for linear elasticity problems, SIAM J. Numer. Anal. 51 (2) (2013) 794–812, http://dx.doi.org/ 10.1137/120874746. [21] P. Wriggers, W.T. Rust, B. Reddy, A virtual element method for contact, Comput. Mech. 58 (6) (2016) 1039–1050. [22] M.L. De Bellis, P. Wriggers, B. Hudobivnik, G. Zavarise, Virtual element formulation for isotropic damage, Finite Elem. Anal. Des. 144 (2018) 38–48. [23] F. Aldakheel, B. Hudobivnik, A. Hussein, P. Wriggers, Phase-field modeling of brittle fracture using an efficient virtual element scheme, Comput. Methods Appl. Mech. Engrg. 341 (2018) 443–466. [24] B. Hudobivnik, F. Aldakheel, P. Wriggers, A low order 3D virtual element formulation for finite elasto–plastic deformations, Comput. Mech. 63 (2) (2019) 253–269. [25] M. De Bellis, P. Wriggers, B. Hudobivnik, Serendipity virtual element formulation for nonlinear elasticity, Comput. Struct. 223 (2019) 106094. [26] P. Wriggers, M. De Bellis, B. Hudobivnik, A Taylor–Hood type virtual element formulations for large incompressible strains, Comput. Methods Appl. Mech. Engrg. 385 (2021) 114021. Computer Methods in Applied Mechanics and Engineering 432 (2024) 117383 20 C. Gatta et al. [27] A. D’Altri, S. de Miranda, L. Patruno, E. Sacco, An enhanced VEM formulation for plane elasticity, Comput. Methods Appl. Mech. Engrg. 376 (2021) http://dx.doi.org/10.1016/j.cma.2020.113663. [28] A. Lamperti, M. Cremonesi, U. Perego, A. Russo, C. Lovadina, A Hu–Washizu variational approach to self-stabilized virtual elements: 2D linear elastostatics, Comput. Mech. 71 (5) (2023) 935–955. [29] M. Nale, C. Gatta, D. Addessi, E. Benvenuti, E. Sacco, An enhanced corotational virtual element method for large displacements in plane elasticity, Comput. Mech. (2024) http://dx.doi.org/10.1007/s00466-023-02437-1. [30] L.L. Yaw, A co-rotational virtual element method for 2D elasticity and plasticity, Internat. J. Numer. Methods Engrg. 125 (6) (2024) e7404. [31] M. Marino, B. Hudobivnik, P. Wriggers, Computational homogenization of polycrystalline materials with the Virtual Element Method, Comput. Methods Appl. Mech. Engrg. 355 (2019) 349–372, http://dx.doi.org/10.1016/j.cma.2019.06.004. [32] M. Pingaro, M.L. De Bellis, E. Reccia, P. Trovalusci, T. Sadowski, Fast statistical homogenization procedure for estimation of effective properties of Ceramic Matrix Composites (CMC) with random microstructure, Compos. Struct. 304 (2023) http://dx.doi.org/10.1016/j.compstruct.2022.116265. [33] N. Sukumar, A. Tabarraei, Conforming polygonal finite elements, Internat. J. Numer. Methods Engrg. 61 (12) (2004) 2045–2066, http://dx.doi.org/10. 1002/nme.1141, Cited by: 353. [34] M. Kraus, A. Rajagopal, P. Steinmann, Investigations on the polygonal finite element method: Constrained adaptive Delaunay tessellation and conformal interpolants, Comput. Struct. 120 (2013) 33–46, http://dx.doi.org/10.1016/j.compstruc.2013.01.017, Cited by: 18. [35] J. Bishop, A displacement-based finite element formulation for general polyhedra using harmonic shape functions, Internat. J. Numer. Methods Engrg. 97 (1) (2014) 1–31, http://dx.doi.org/10.1002/nme.4562, Cited by: 93. [36] G. Manzini, A. Russo, N. Sukumar, New perspectives on polygonal and polyhedral finite element methods, Math. Models Methods Appl. Sci. 24 (8) (2014) 1665–1699, http://dx.doi.org/10.1142/S0218202514400065, Cited by: 115. [37] A. Francis, A. Ortiz-Bernardin, S.P. Bordas, S. Natarajan, Linear smoothed polygonal and polyhedral finite elements, Internat. J. Numer. Methods Engrg. 109 (9) (2017) 1263–1288, http://dx.doi.org/10.1002/nme.5324, Cited by: 64. [38] S. Ghosh, R. Mallett, Voronoi cell finite elements, Comput. Struct. 50 (1) (1994) 33–46, http://dx.doi.org/10.1016/0045-7949(94)90435-9. [39] S. Ghosh, K. Lee, S. Moorthy, Multiple scale analysis of heterogeneous elastic structures using homogenization theory and voronoi cell finite element method, Int. J. Solids Struct. 32 (1) (1995) 27–62, http://dx.doi.org/10.1016/0020-7683(94)00097-G, Cited by: 344. [40] S. Ghosh, K. Lee, P. Raghavan, A multi-level computational model for multi-scale damage analysis in composite and porous materials, Int. J. Solids Struct. 38 (14) (2001) 2335–2385, http://dx.doi.org/10.1016/S0020-7683(00)00167-0, Cited by: 309. [41] M. Groeber, S. Ghosh, M.D. Uchic, D.M. Dimiduk, A framework for automated analysis and simulation of 3D polycrystalline microstructures. Part 1: Statistical characterization, Acta Mater. 56 (6) (2008) 1257–1273, http://dx.doi.org/10.1016/j.actamat.2007.11.041, Cited by: 241. [42] P.L. Bishay, S.N. Atluri, Trefftz-Lekhnitskii Grains (TLGs) for efficient direct numerical simulation (DNS) of the micro/meso mechanics of porous piezoelectric materials, Comput. Mater. Sci. 83 (2014) 235–249, http://dx.doi.org/10.1016/j.commatsci.2013.10.038. [43] J. Reinoso, M. Paggi, A consistent interface element formulation for geometrical and material nonlinearities, Comput. Mech. 54 (2014) 1569–1581. [44] D. Addessi, P. Di Re, E. Sacco, Micromechanical and multiscale computational modeling for stability analysis of masonry elements, Eng. Struct. 211 (2020) 110428. [45] G. Borino, F. Parrinello, A symmetric tangent stiffness approach to cohesive mechanical interfaces in large displacements, Int. J. Comput. Methods Eng. Sci. Mech. 23 (6) (2022) 551–566. [46] E. Sacco, Revisiting Some Modeling of Damage-Friction-Dilatancy Coupling in Cohesive Interfaces, Available at SSRN 4664241. [47] S. Marfia, E. Monaldo, E. Sacco, Cohesive fracture evolution within virtual element method, Eng. Fract. Mech. 269 (2022) 108464. [48] M.F. Benedetto, A. Caggiano, G. Etse, Virtual elements and zero thickness interface-based approach for fracture analysis of heterogeneous materials, Comput. Methods Appl. Mech. Engrg. 338 (2018) 41–67. [49] S. Marfia, E. Monaldo, E. Sacco, Cohesive fracture evolution within virtual element method, Eng. Fract. Mech. 269 (2022) http://dx.doi.org/10.1016/j. engfracmech.2022.108464. [50] M. Bruggi, P. Venini, Modeling cohesive crack growth via a truly-mixed formulation, Comput. Methods Appl. Mech. Engrg. 198 (47–48) (2009) 3836–3851, http://dx.doi.org/10.1016/j.cma.2009.08.018. [51] U. De Maio, F. Greco, L. Leonetti, P. Nevone Blasi, A. Pranno, A cohesive fracture model for predicting crack spacing and crack width in reinforced concrete structures, Eng. Fail. Anal. 139 (2022) http://dx.doi.org/10.1016/j.engfailanal.2022.106452. [52] M. Ostoja-Starzewski, Microstructural Randomness and Scaling in Mechanics of Materials, CRC Press, Taylor & Francis Group, 2007. [53] M. Ostoja-Starzewski, Macrohomogeneity condition in dynamics of micropolar media, Arch. Appl. Mech. 81 (7) (2011) 899–906, http://dx.doi.org/10. 1007/s00419-010-0456-1. [54] M. Ostoja-Starzewski, Material spatial randomness: From statistical to representative volume element, Probab. Eng. Mech. 21 (2) (2006) 112–132, http://dx.doi.org/10.1016/j.probengmech.2005.07.007. [55] Z. Khisaeva, M. Ostoja-Starzewski, On the size of RVE in finite elasticity of random composites, J. Elasticity 85 (2) (2006) 153–173, http://dx.doi.org/ 10.1007/s10659-006-9076-y. [56] M. Pingaro, E. Reccia, P. Trovalusci, R. Masiani, Fast statistical homogenization procedure (FSHP) for particle random composites using virtual element method, Comput. Mech. 64 (1) (2019) 197–210, http://dx.doi.org/10.1007/s00466-018-1665-7. [57] M. Pingaro, E. Reccia, P. Trovalusci, Homogenization of random porous materials with low-order virtual elements, ASCE-ASME J. Risk Uncertain. Eng. Syst. B: Mech. Eng. 5 (3) (2019) http://dx.doi.org/10.1115/1.4043475. [58] V. Eremeyev, J.-F. Ganghoffer, V. Konopińska-Zmysł owska, N. Uglov, Flexoelectricity and apparent piezoelectricity of a pantographic micro-bar, Internat. J. Engrg. Sci. 149 (2020) http://dx.doi.org/10.1016/j.ijengsci.2020.103213. [59] N. Mawassy, H. Reda, J.-F. Ganghoffer, V. Eremeyev, H. Lakiss, A variational approach of homogenization of piezoelectric composites towards piezoelectric and flexoelectric effective media, Internat. J. Engrg. Sci. 158 (2021) http://dx.doi.org/10.1016/j.ijengsci.2020.103410. [60] P. Trovalusci, M.L. De Bellis, M. Ostoja-Starzewski, A. Murrali, Particulate random composites homogenized as micropolar materials, Meccanica 49 (11) (2014) 2719–2727. [61] P. Trovalusci, M. Ostoja-Starzewski, M. De Bellis, A. Murrali, Scale-dependent homogenization of random composites as micropolar continua, Eur. J. Mech. A Solids 49 (2015) 396–407, http://dx.doi.org/10.1016/j.euromechsol.2014.08.010. [62] E. Reccia, M.L. De Bellis, P. Trovalusci, R. Masiani, Sensitivity to material contrast in homogenization of random particle composites as micropolar continua, Composites B 136 (2018) 39–45. [63] R. Smit, W. Brekelmans, H. Meijer, Prediction of the mechanical behavior of nonlinear heterogeneous systems by multi-level finite element modeling, Comput. Methods Appl. Mech. Engrg. 155 (1–2) (1998) 181–192, http://dx.doi.org/10.1016/S0045-7825(97)00139-4. [64] C. Miehe, J. Schröder, J. Schotte, Computational homogenization analysis in finite plasticity simulation of texture development in polycrystalline materials, Comput. Methods Appl. Mech. Engrg. 171 (3–4) (1999) 387–418, http://dx.doi.org/10.1016/S0045-7825(98)00218-7. [65] R. Hill, Elastic properties of reinforced solids: Some theoretical principles, J. Mech. Phys. Solids 11 (5) (1963) 357–372. [66] M. Sena, M. Ostoja-Starzewski, L. Costa, Stiffness tensor random fields through upscaling of planar random materials, Probab. Eng. Mech. 34 (2013) 131–156, http://dx.doi.org/10.1016/j.probengmech.2013.08.008. [67] X. Du, M. Ostoja-Starzewski, On the size of representative volume element for Darcy law in random media, Proc. R. Soc. A 462 (2074) (2006) 2949–2963, http://dx.doi.org/10.1098/rspa.2006.1704. Computer Methods in Applied Mechanics and Engineering 432 (2024) 117383 21 C. Gatta et al. [68] C. Talischi, G. Paulino, A. Pereira, I. Menezes, PolyMesher: A general-purpose mesh generator for polygonal elements written in Matlab, Struct. Multidiscip. Optim. 45 (3) (2012) 309–328, http://dx.doi.org/10.1007/s00158-011-0706-z. [69] M. Pingaro, M.L. De Bellis, P. Trovalusci, R. Masiani, Statistical homogenization of polycrystal composite materials with thin interfaces using virtual element method, Compos. Struct. 264 (2021) 113741. [70] L. Beirão Da Veiga, F. Brezzi, L.D. Marini, A. Russo, The virtual element method, Acta Numer. 32 (2023) 123–202, http://dx.doi.org/10.1017/ S0962492922000095. [71] P. Wriggers, B. Hudobivnik, F. Aldakheel, A virtual element formulation for general element shapes, Comput. Mech. 66 (4) (2020) 963–977, http: //dx.doi.org/10.1007/s00466-020-01891-5. [72] E. Artioli, L. Beirão da Veiga, C. Lovadina, E. Sacco, Arbitrary order 2D virtual elements for polygonal meshes: part I, elastic problem, Comput. Mech. 60 (3) (2017) 355–377, http://dx.doi.org/10.1007/s00466-017-1404-5. [73] E. Sacco, J. Toti, Interface elements for the analysis of masonry structures, Int. J. Comput. Methods Eng. Sci. Mech. 11 (6) (2010) 354–373. [74] D. Addessi, C. Gatta, S. Marfia, E. Sacco, Multiscale analysis of in-plane masonry walls accounting for degradation and frictional effects, Int. J. Multiscale Comput. Eng. 18 (2) (2020). [75] D. Addessi, P. Di Re, C. Gatta, E. Sacco, Non-uniform TFA reduced multiscale procedure for shell-3D modeling of periodic masonry structures, Mech. Res. Commun. 130 (2023) 104122. [76] G. Alfano, E. Sacco, Combining interface damage and friction in a cohesive-zone model, Internat. J. Numer. Methods Engrg. 68 (5) (2006) 542–582. [77] M. Paggi, P. Wriggers, A nonlocal cohesive zone model for finite thickness interfaces–Part II: FE implementation and application to polycrystalline materials, Comput. Mater. Sci. 50 (5) (2011) 1634–1643. [78] Z. Gong, W. Zhao, K. Guan, P. Rao, Q. Zeng, J. Liu, Z. Feng, Influence of grain boundary and grain size on the mechanical properties of polycrystalline ceramics: Grain-scale simulations, J. Am. Ceram. Soc. 103 (10) (2020) 5900–5913. [79] Y. Zhao, Q. Song, H. Ji, W. Cai, Z. Liu, Y. Cai, Multi-scale modeling method for polycrystalline materials considering grain boundary misorientation angle, Mater. Des. 221 (2022) 110998. [80] P. Trovalusci, M.L. De Bellis, R. Masiani, A multiscale description of particle composites: From lattice microstructures to micropolar continua, Composites B 128 (2017) 164–173. Computer Methods in Applied Mechanics and Engineering 432 (2024) 117383 22