Homogenization of carbon/polymer composites with anisotropic distribution of particles and stochastic interface defects

The main objective is to investigate an effect of anisotropic distribution of the reinforcing particles in a cubic representative volume element (RVE) of the carbon–polymer composite including stochastic interphases on its homogenized elastic characteristics. This is done using a probabilistic homogenization technique implemented using a triple approach based on the stochastic perturbation method, Monte Carlo simulation as well as on the semi-analytical approach. On the other hand, the finite element method solution to the uniform deformations of this RVE is carried out in the system ABAQUS. This composite model consists of two neighboring scales–the micro-contact scale relevant to the imperfect interface and the micro-scale—having 27 particles inside a cubic volume of the polymeric matrix. Stochastic interface defects in the form of semi-spheres with Gaussian radius are replaced with the interphase having probabilistically averaged elastic properties, and then such a three-component composite is subjected to computational homogenization on the microscale. The computational experiments described here include FEM error analysis, sensitivity assessment, deterministic results as well as the basic probabilistic moments and coefficients (expectations, deviations, skewness and kurtosis) of all the components of the effective elasticity tensor. They also include quantification of anisotropy of this stiffness tensor using the Zener, Chung–Buessem and the universal anisotropy indexes. A new tensor anisotropy index is proposed that quantifies anisotropy on the basis of all not null tensor coefficients and remains effective also for tensors other than cubic (orthotropic, triclinic and also monoclinic). Some comparison with previous analyses concerning the isotropic case is also included to demonstrate the anisotropy effect as well as the numerical effort to study randomness in composites with anisotropic distribution of reinforcements and inclusions.


Introduction
It is widely known that the interface defects appearing between different phases of composites (like a matrix and its reinforcement) play a crucial role in reliability, durability and failure of these materials [1][2][3]. Generally, structural defects in composites may reflect some manufacturing imperfections [4], following remarkable residual thermal stresses for instance, and may exhibit spherical or ellipsoidal shapes, at least for polymeric materials [5], resulting from the well-documented matrix cavitation processes [6]. Therefore, it seems reasonable that interface inaccuracies or imperfections can be approximated by semicircular or semi-elliptical shapes in 2D models of the composites (or by semi-spherical or semi-ellipsoidal shapes in 3D analysis) [7,8]. Particulate or fibrous reinforcements usually do not contain such imperfections unlike the matrices, whereas the largest thermal stresses appear at the reinforcement-matrix boundary, so that the largest concentration of these imperfections is expected at the interface and they are geometrically located in the matrix region. Different numerical models are available in the literature to include this effect in numerical simulation of the composite-(i) from contact finite elements [9], (ii) through the interface elements [10] or the set of disconnecting spring elements up to (iii) the interphase [8,[11][12][13] being a thin film in-between traditional composite constituents modeled with 2 or 3D solid elements having material characteristics remarkably weaker than these defined for the matrix. One needs to notice a variety of analytical estimates of the effective elasticity tensor components for the composites including some interface imperfections [14][15][16]. Computational cost of numerical simulation of such a composite remarkably increases when compared to ideal two-component structure, especially when more than a single reinforcement is inserted into the RVE and when uncertainty analysis is to be carried out. Therefore, a homogenization method [13,14] is recommended to reduce this complexity by scale reduction in complex composite structures from a three-scale approach to traditional micro-macroanalysis where two separate scales are really necessary. Such a reduction gains specific importance when anisotropic packaging of the particles or fibers is assumed into the RVE [17][18][19][20], because inserting additional imperfect interfaces or interphases dramatically complicates meshing of the entire structure; it is also really necessary when statistical simulation techniques [21,22] are employed together with the complex finite element method (FEM) composite model.
A complete lack of information concerning effective behavior of the particle-reinforced composite with anisotropic distribution of the particles having randomly imperfect contact with the surrounding matrix was the main motivation to complete the study presented below. Three concurrent probabilistic methods have been used-Monte Carlo simulation [21], iterative stochastic perturbation technique (called here ISFEM after corresponding FEM implementation) [23] as well as the semi-analytical method [24,25], all based on the finite element method experiments to verify an influence of the interphase with random volume fractions of the defects on the effective characteristics of such a composite. We apply for this purpose Gaussian parameter w with the given expectation E[w] and standard deviation σ (w) to compute the first four probabilistic moments and coefficients of the effective tensor and to verify whether it also can have Gaussian distribution or not. A choice of this parameter originates with experimentally driven uncertainty in the micro-cavities appearing in the matrix and localized at the particle-matrix interface. They are modeled here as semi-spheres with Gaussian radius, whereas the fiber-reinforced structure models use semicircular or semi-elliptical approximations [7,8]. A homogenization method is applied here either to a composite with isotropic distribution of the particles and also for accidental anisotropic location of carbon black particles to verify an influence of this variation. A decisive part of this study is the weighted least squares method (WLSM) approximation of the homogenized tensor components as polynomial functions of this w of numerically optimized order. Three criteria are employed here-maximization of correlation of the FEM experiments series to the designed polynomial as well as minimization of the variance and error in this fitting procedure. Let us note that mesh generation for the composite with anisotropically distributed particles was the very challenging problem itself having even some literature sources [26]. Determination of the optimal polynomial representations of the effective tensor components versus parameter w enables for an analytical derivation of these tensor probabilistic moments as the functions of E[w] and σ (w), which is also contained in this work. Sensitivity analysis of the homogenized tensor and the FEM error analysis have been carried out and discussed in detail in parallel to the main probabilistic analysis, which can further lead to optimal design of the new composites [27,28]. The largest potential of this work is in application of nonlinear constitutive models [29] of the polymer especially, as well as statistical experimental information on the interfaces [30], which can give the input to further probabilistic studies. The entire study is a next milestone to carry out reliability assessment of various particle-reinforced materials and structures with the use of a homogenization method, without a need of very precise multi-scale meshing of the composite structural elements.  [13] 2 Mathematical model Let us consider a particle-reinforced composite in 3D Euclidean space indexed with Cartesian coordinates and denote it by Y . Its representative volume element (RVE) is denoted by and consists of three disjoint subsets corresponding to spherical particles ( p ), to the matrix ( m ) as well as to the interface defects ( d ) such that (1) A perfect and load-independent contact is assumed in between all these subsets during further deformation process. Let ∂ denote the outer surface of the RVE and let p ∩ ∂ = ∅. Without any loss of generality, the RVE under consideration contains r separate reinforcing particles non-uniformly distributed in the volume of having all the same distribution of the defects. It is assumed that the particles and the matrix work in linear elastic reversible regime, where Young's moduli and Poisson ratios are defined as The interface defects subset can be further decomposed into a series of disjoint single semi-bubbles in the following way: We assume further that a single interface defect d(i, j) has semi-spherical shape, is localized on the given particle surface and is directed outwards from the particle center (Fig. 1). A subscript i stands for the defects number on the given particle-matrix interface (1 ≤ i ≤ n), while j = 1, ..., r indices the particle number. Further, we assume that interface defects throughout the given particle-matrix interface (their radii R (i, j) ) are statistically scattered according to the Gaussian distribution defined by its first two moments E R (i, j) and Var R (i, j) .
Multi-scale homogenization method presented in the next section is based on two separate steps: (i) reconstruction of the interphase (third component) as the new artificial material fully separating particles from the surrounding matrix, (ii) probabilistic averaging of the defects inside the interphase where they belong to (micro-contact scale), (iii) probabilistic homogenization of the entire three-component composite (microscale). This interphase is a thin layer around any particle, whose constant thickness equals to an upper bound on statistical population of the defects radii according to the well-known three-sigma rule: Finally, we propose a probabilistic averaging method to recover effective elastic properties of the interphase reduced in addition to the original properties of the matrix. It starts from deterministic formula describing Young's modulus of the given interphase as volumetric average of the defects and of the matrix belonging to its volume i . There holds Recalling that E d = 0 and substituting volume ratio of the defects inside the interphase as the parameter w one can notice that the first two moments of the variable E i are equal to where z indices three constituents (interphase i, particle p and matrix m). Calculation of the effective elasticity tensor components proceeds by using of the following formula: which represents deformation energies of the real (at the left-hand side) and also of the homogenized (on the right-hand side) composite RVE. The R.H.S. in this equation can be replaced with more detailed deformation energy distribution into the particle, the interphase and the matrix as Determination of the first four probabilistic moments of the effective elasticity tensor proceeds directly from the their definitions and returns consecutively It is possible to scale the deformations ε αβ to make the L.H.S. integral equal to 1 and then, the moments in Eqs. (11)(12)(13) just equal to the moments of deformation energy of the real composite. On the other hand, it is well known that elasticity tensor for the anisotropic continuum has definitely more complex form than in Eq. (7) and one writes it in a form of the fourth-or second-rank tensors as The effective elasticity tensor for the homogenized composite, which is fully anisotropic (after non-uniform distribution of the isotropic particles into the isotropic matrix plus interphases) is calculated here by a substitution of the results of the strain energy U obtained from the FEM simulations of 21 independent unit displacements of the RVE. Three of these are extensional (u 1 , u 2 , u 3 = 1), also 3-transverse (u 12 , u 13 , u 23 = 1) and 15 remaining constitute coupling displacements u i j u kl = 1 for i, j, k, l ∈ 1, 2, 3 , i = j = k = l and further inserted into Eq. (8), e.g., taking boundary as u 1 = 1 and u i j = 0 for i, j = 1 one can obtain ε 11 = 1 where all the remaining strains simply vanish.
Finally, we find numerically polynomial representations of the effective elasticity tensor components (and also of the deformation energy stored in the RVE), e.g., and their Least Squares Method fittings of the minimum order are equivalent to n=4 (5 and 6 also), so that analytical derivation of the first two probabilistic moments of this tensor in case of Gaussian uncertainty in w (with mean the value E[w] and standard deviation σ (w)) yields and also (with no summation on α, β, χ , δ = 1, 2, 3 at the R.H.S.) These equations have been derived assuming that each particle is surrounded by the interphase of the same thickness and the same elastic properties reduction parameter w. This thickness does not appear as an explicit additional parameter here but affects computational FEM experiments leading to the specific LSM approximation of the effective elasticity tensor. It is possible to insert different Gaussian random variables w i here, whose basic probabilistic moments may be different for various particles interphases but then a multi-parametric LSM fitting is necessary and additionally, the cross-correlations of the effective tensors are to be calculated also. Numerical experiments provided here show that both isotropic and anisotropic distributions of the particles result in the same order polynomial representations of expected values and variances; some differences are obtained for the coefficients a αβγ δ(i) only. On the other hand, when the interphase thickness approaches 0, then the expected values of effective tensor components become constant and equal to a αβγδ(0) ; the expectation starts to be independent of the parameter w in the limiting case. The variances of the homogenized tensor components remarkably all tend to 0, so that these components themselves become simply deterministic. A degree of the LSM fitting polynomial may fluctuate when the interphase thickness tends to 0; however, this cannot change the observations collected above. Let us also note that higher-order central probabilistic moments and 5th-or 6thorder polynomial expressions and their moments have been omitted for a brevity of the presentation. Therefore, one can demonstrate analytically following considerations attached in the Appendix, that effective elasticity tensor components cannot have Gaussian distribution while randomizing interface defects ratio according to this distribution. Therefore, the first two moments characterization of this tensor does not identify it uniquely and one needs higher-order stochastic computational technique to reliably determine its constituents.

Deterministic computational analysis
Computational experiments with homogenization of the particle-reinforced composite with anisotropic distribution of the particles inside its RVE and imperfect interfaces consist of the following tasks: 1. computational deformations of the RVE under uniaxial and biaxial extension corresponding to the randomized mechanical properties of the interphase E i , ν i , 2. relative computational error relevant to the uniaxial and biaxial extension of the RVE corresponding to the boundary conditions imposed on the boundary problem, 3. sensitivity coefficients S C (eff) αβχ δ with respect to the consecutive volume fractions of the interface defects w, 4. optimal order determination for the polynomial response functions of the effective tensor components C (eff) αβχ δ = C (eff) αβχ δ (w), 5. final calculation of the expected values, coefficients of variation, skewness and kurtosis of C (eff) αβχ δ . Computer model of the composite represents a homogenization method of three different representative volume elements (RVEs) with three phases and 27 isotropic (3rd model) or anisotropic (1st and 2nd model) distributions of the spherical particles surrounded by the same spherical interphases (Fig. 2). Such a distribution of the reinforcing particles results in an isotropic (1st model) or anisotropic (2nd and 3rd models) effective elasticity tensor. RVEs have dimensions of [3δ, 3δ, 3δ], where δ is unit interval and they are composed of a polymeric matrix, a spherical carbon black C 60 particle and an interphase of constant thickness in between these two materials. Particle location is chosen in an accidental way and numerical algorithms presented in [31] or [32] are omitted here principally because this study focuses on an anisotropic representation of the homogenized composite and its impact on the effective tensor. A number of the particles representing an internal structure of the composite has been chosen as 27 because (1) it ensures the best balance between quality of its internal composition and mesh quality limited by the computational resources, (2) it allows adequate representation of internal structure in cubic composite and (3) it is very close to the number of monodisperse particle chosen in other studies [32]. These particles are not restricted inside the RVE with an exception to the 3rd model, whose internal composition ensures a perfect regular cubic distribution. Initial mechanical properties of all phases are adopted as follows, for the matrix asE m = 4.0 MPa, ν m = 0.34, for the CB particle as E p = 10.0 GPa, ν p = 0.3, and vanishing both for the interface defects. All computational experiments with the RVE have been carried out in the FEM system ABAQUS and consist of a series of the cell problems solved for deformation of 21 different uni-and bidirectional stretches as well as shear deformations. A set of all the required boundary conditions are shown in Fig. 3, where they are grouped into six different cases, namely three uniaxial tensions (1) and shear deformations (2), three biaxial tensions (3) and shear deformations (4), tri-axial tension and corresponding shear deformations (5) a complete set of all six shear deformations and extensions (6). Colors used in this graph represent a magnitude of displacement, whose maximal value reaches 6.0 and is obtained from an 100% uniaxial stretch coupled with 100% shear in the same direction. A limitation or averaging of the amount of unit displacements is not applied here because a result of homogenization is a fully anisotropic stiffness tensor; it is compared for three RVEs signifying structuring, anisotropy or clustering of the particles. Each of the RVEs has been discretized separately with hexahedral 8-noded C3D8 1st-order finite elements, hexahedral 20-noded C3D20, and also using tetrahedral 10-noded C3D10 2nd-order finite elements. The 1st model is composed of hexahedral 1st-order elements with amount of 239.488 (1st model), and the 3rd model-36.504 2nd-order hexahedral elements, while the 2nd model is combined with 9.504 hexahedral and 174.044 tetrahedral 2nd-order elements (2nd model). The type and amount of finite elements in each model is selected to obtain the best possible numerical accuracy with acceptable computational effort, which is verified using FEM error analysis presented and discussed in Figs. 4 and 5. The amount of finite elements (T F E ) is not constant for these models, because non-uniform distribution of the particles in the 2nd model needs tetrahedral mesh of the matrix. Although anisotropic, the first model prevents some structuring of particles and, therefore, hexahedral mesh for matrix is possible. Regardless of the finite element type and model number, a tie contact must be placed between the interphase and the matrix to ensure deformation continuity. Uniform extensions and shear deformations in all models are applied to use energybased homogenization method (see Fig. 2) to calculate effective elasticity tensor components and, finally, to contrast its probabilistic characteristics resulting from three independent numerical techniques to verify the influence of (1) particles anisotropy and (2) uncertainty in interface defects. The first two models are composed of non-homogeneously distributed particles, and their elasticity tensor has 21 independent coefficients. These coefficients can be further reduced to 6 representative groups indicated on the bottom left part of the Fig. 2 as they have been computed for very similar boundary conditions (see Fig. 3). The first model, additionally to this geometrical anisotropy, has also two sets of clustering particles with an amount of four and two. Situation is totally different for the third model with a uniform distribution of the reinforcing particles, whose deformation can be entirely represented by only three independent elasticity tensor components (C 11 , C 44 and C 12 ).
Further computations have been carried out entirely using symbolic algebra system MAPLE and are focused on determination of the probabilistic parameters of the RVE effective tensor components C (eff) with respect to the volume fraction of interface defects in the interphase using polynomial bases. These contain determination and optimization of the polynomial response function of the effective tensor components using the WLSM with Dirac-type distribution of weights [1,1,1,1,1,10,1,1,1,1,1]. Further, they concern computation of the basic probabilistic characteristics of this tensor with input Gaussian uncertainty in the volume fraction of interface defects in the interphase w. FEM simulations for each mean value of w are performed for its ±10% neighborhood, whereas the Weighted Least Squares Method (WLSM) is implemented in its nonlinear version, where additional parameters are the WLSM increments number and polynomial order, while statistical optimization is done via correlation maximization as well as variance and RMS error minimization of the WLSM-a chosen number of the increments equals 15 and polynomial order equals 4, 5 or 6 depending on the effective tensor component. Interestingly, the WLSM for polynomials above the 6th order returns much higher RMS error and variance with only slightly decreased correlation when comparing them with the given optimal order. Stochastic simulations are performed separately for three RVEs and for an increasing expected value of the volume fraction of interface defects w by using the following three independent methods: semi-analytical, stochastic perturbation of second and fourth order as well as Monte Carlo simulation with 250.000 random trials. This enables a triple comparison of (a) statistical scattering in anisotropic and cubic effective tensor, (b) the efficiency of the applied probabilistic methods and their coincidence, as well as (c) a verification how an increase in volume fraction of the interface defects affects uncertainty in the resulting effective tensor. These computations include studies of the expected values ( 6 representative tensor components C (eff) fully illustrating its uncertain response to the random volume fraction of interface defects w. It should be additionally noted that isotropic effective tensor of the third model has only three components C 22 , C 55 and C 12 .
Since the major focus of this paper is comparison of two anisotropic and isotropic random effective tensor, all figures (Figs. 4-45) are composed in the same way, where the given symbols represent the model type, i.e., asterisk-first model (anisotropic), circle-second model (anisotropic), rectangle-third model (cubic) and diamond-analytical result (cubic). The colors have a meaning specific to each set of figures, e.g., for effective tensor coefficients (Figs. 6,7,8,9,10,11) they distinguish between the coefficients or for all the probabilistic results on the graph left side (Figs. 22a-45a); they also allow to distinguish the expected values of the volume fraction of interface defects E (w).  Numerical analysis of the FEM error of the total strain energy for all the RVE models, which enters computer analysis, is here twofold. First of all, computational FEM error obtained for different types of finite elements is compared for uniaxial tension with respect to the total amount of the FEs in the RVE denoted by T F E .This allows to determine the most suitable element type for each model (Fig. 4). Secondly, an error for the chosen element types is further checked for all the representative groups of components (distinguished by colors on the graph) in relation to the total number of variables in the FEM solver T var (Fig. 5). This is done (a) to investigate an efficiency of the particular finite element type for all the necessary sets of boundary conditions in the FEM homogenization tests and (b) to enable a proper choice of a number of these elements (variables) ensuring satisfactory convergence of the results, separately for each model. Thanks to these two analyses, it is possible to contrast results of all the three models eliminating an impact of discretization intensity and element type on the results of the strain energies U , which are then converted into the effective tensor components C (eff) . A relative FEM error for both figures has been calculated in addition to the computation of the highest element amount for each element type according to the following formulas: is the total strain energy for the ith FEM simulation, U T F E max stands for the total strain energy for the FEM simulation with a highest amount of FEs of the specific type, and these are replaced by the U T var(i) and U T var max in Fig. 5, where T F E is substituted by T var .
Computations contained in Fig. 4 include finite elements of (a) hexahedral first order with full integration C3D8 (hex 1st), (b) hexahedral second order with full integration C3D20 (hex 2nd), (c) tetrahedral first-order C3D4 (tet 1st), (d) tetrahedral second-order C3D10 (tet 2nd) or (e) and (f) a combination of these two types of elements of the same order (hex&tet 1st and hex&tet 2nd). This graph shows firstly a very small error of the FEM simulations for all the considered element types and amount almost never exceeding 2% except tetrahedral elements of the first order and generally less than 12 for the T F E adopted in further computations (Fig. 5). The second-order formulation of the FEs ensures a slightly better performance, which is especially remarkable in the third model, where the error even does not reach 2 . This performance is not so high for the anisotropic models, but still at least 30% better than for the first-order elements. An additional observation is that the choice of tetrahedral elements results in a higher relative error (as they have slower convergence) and that they are not recommended for the domains exhibiting axisymmetric geometry (such as the spherical carbon black particle or an interphase), since they result in the error close to 18%. A major reason for this is that they are, in contrast to the hexahedral elements, implemented by unstructured mesh generators. This is why for both the first and the third model, the hexahedral elements are chosen with of almost 240.000 1st-order elements (for the first model) and 35.000 2nd-order elements (for the second model).
On the other hand, it should be noticed that in order to enable a hexahedral mesh, a part (e.g., the matrix) has to be in many cases partitioned to gain a structured regions, where automatic generators serve efficiently. This is the case of both the first and the third models (both partitioned by a couple of dozen planes) and implies that for highly complex geometries, especially in a 3D stress state, the model must be cut into the few smaller parts. This, in turn, (1) restricts FE usage to the very small elements, (2) considerably increases the modeling and computation time and (3) is difficult or unavailable for an automatic generation. This may prevent usage of the 2nd-order hex formulation having higher restrictions on geometry, as in case of the 1st model and where manual partitioning may be even not possible (for totally unstructured parts). Such a situation takes place for the second model, whose particles have absolutely no structuring and, therefore, the hexahedral mesh cannot be applied for its matrix. It is, however, still better to preserve the hexahedral mesh for the parts with axisymmetric geometry (particle and the interphase), while assigning a tetrahedral mesh to the matrix. This is well justified by Fig. 3, where such an approach gives almost the results about two times more precise than the tetrahedral mesh for the 1st and the 2nd-order elements, particularly for the small T F E . On the basis of the above, a 2nd-order hexahedral mesh of amount of almost 10.000 is assigned in this model for the particle and the interphase, which is supplemented by the 2nd-order tetrahedral mesh of more than 170.000 elements for the matrix; similarly to the remaining two models, such a choice ensures a relative error much below 1 . From Fig. 4, it can be also concluded that FEM numerical error sharply decreases together with an increase in T F E and that the 3rd and the least complex model has the smallest error for all the element types. This asset, however, vanishes together with an increasing T F E , where the error asymptotically goes to 0 (consistently with the FEM theory). Therefore, all the models are convergent, but the anisotropic ones require more than 10 times higher T F E to preserve the same error level like isotropic composite.
Further, computational time and physical memory usage of the FEM solution are also of paramount importance and they are usually several or even more than ten times higher for the second-order formulation for the same T F E . This difference is much smaller when contrasting these formulations by the total amount of variables in the FEM T var , and this relation is shown in Fig. 5. Nonetheless, it should be still considered that even when plotting the error in relation to the T var , tetrahedral elements have from 3 up to 4 times longer preprocessing than the hexahedral mesh, while the 2nd-order elements consume 20-50% more processor time and the 2nd-order formulation use up to 80% more physical memory for the same T var (not 10 times or more as for the T F E ). Figure 5 firstly shows that even for 10.000 finite elements of the given type an overall computational error is very small and smaller or equal than 4 . This error is very similar for each model for all the given deformations except the first model, where the shear deformations have up to three times smaller error for a low T var . Secondly, the chosen element types are also highly convergent to a 0 error for all but the second model and vanish almost completely for the optimum T var representing the mesh density used in further computations. This second model is here fluctuating a little bit together with an increase in T var , which may be caused by the use of an unstructured tetrahedral mesh for the matrix. Such a mesh is possibly to be generated in a different way or from a different point even for small changes in mesh seeds. Still, these variations are in the range of several and would not affect the computation efficiency of the FEM homogenization problem solution. The last remark is that 2nd-order hexahedral elements in the third model give a higher amount of variables than for the anisotropic ones. This firstly results from the rich formulation of such an element and partitioning of the RVE cell as well as rendering solver problems for small amount of elements. It additionally has a connection with a method of defining RVE geometry for the FEM solver, a little different for this model than for the anisotropic ones due to its internal symmetries-in this case the RVE is composed of 27 equal cells.
The RVEs deformation energies for different models necessary in homogenization method are converted first to the 21 effective tensor components, then grouped into the 6 representative sets and finally plotted in relation with w. This is done separately for five mean values of the volume fraction of interface defects in the interphase having a magnitude equal to (a) 0.10, (b) 0.20, (c) 0.30, (d) 0.40, and (d) 0.50 in their 10% neighborhood giving 55 separate FEM computations of each deformation into a single FEM model, 1.155 for each model and 3.465 in total. These effective tensor components are depicted for all the effective tensor groups in Figs. 6, 7, 8, 9, 10, 11 and 12 as separate discrete points per each FEM test-they are plotted as the functions of the parameter w for three different RVE models. Additionally, a result of an analytical approach calculated according to Christensen [33] is plotted for C 11 , C 44 and C 12 with w = 0 and marked using a brown diamond to validate further a correctness of the computationally derived effective tensor components. Null w represents a special case when properties of the interphase are equal to the ones of the matrix-the interphase is not existent. Results of the third model are available only for Figs. 5, 6 and 7, where the three independent coefficients C 11 , C 44 and C 12 are repeated for other axes (e.g., C 11 ≡ C 22 ≡ C 33 ), which is done for better transparency of the legend. Such a relation is an analytical derivative of the isotropic effective tensor (its polynomial response representation), appropriate for the third model. The effective tensor components C eff αβ are generally smooth, continuous and decreasing with an increase in w. The ones where the third-model tensor is available (Figs. 5, 6, 7) show quite satisfactory convergence of the analytical results with the computational models for w = 0. This is especially true for the isotropic (third) model, which directly represents the analytical model and gives the closest results for C 44 and C 12 ; it is a little bit more distant then the second model for C 11 , but the difference is smaller than 1% and thus negligible. The anisotropic models, being close to the analytical approach for the diagonal effective components (group of C 11 and C 12 ), exhibit much higher shearing effective coefficients for w = 0. It must be underlined here that these two models may, and even should, result in a slightly different effective components than then analytical ones. This difference from the third (analytical) model represents an influence of anisotropy and clustering. These models are treated here as more accurate in the sense that they better reflect micro-geometry of this composite. A second remark is that the simple regular isotropic (third) model overestimates the tensional effective components especially Effective tensor components corresponding to a pure shear for three spatial distributions of particles in the RVE for larger w, while generally underestimating the shearing effective components (the higher w, the smaller underestimation). Importantly, the anisotropic models are much more affected by the interphase defects than the cubic one, so that their effective tensor components in each group are reduced in the same way by these defects and keep close to each other. This means that anisotropy has rather global influence on the effective tensor, affects especially the shearing components and makes the composite much more compliant to this interphase.
Clustering of particles included in the first model does not have a predicted influence on the effective tensor since it increases the magnitude of this tensor components up to 2% while compared to the second model with anisotropy and no clustering. This difference, being negligible for the composite without the interface defects, increases together with w and is the highest for the component corresponding to the uniaxial tension. This effective tensor is still lower than computed for the isotropic particles location, but such results are not consistent with practical observations, where clustering is one of the main reasons for softening of particle-reinforced composites. Stress distributions and magnitudes resulting in different models (Fig. 8) show almost four times higher intensity for the model with clustering localized of course within the interphase. There is no doubt that the interphase very early enters the nonlinear regime and accumulates damage with such high reduced stress intensity. Therefore, this result is limited to the elastic region and clustering certainly softens such a composite (reduce the magnitudes of effective tensor components) when more accurate constitutive model is chosen like some hyper-elastic stress-strain experimental curve available for the polymer matrix possibly including damage or softening of the interphase. A comprehensive comparison of the stress intensity appearing in the models is postponed here for the presentation brevity.
It looks that this composite is the stiffest in uniaxial tension, where a magnitude of the effective tensor components reaches up to 6.9 · 10 6 Pa for w = 0.0 and is reduced by the defects with w = 0.5 to 5.6 · 10 6 Pa (Fig. 6). Its resistance to a single-axial shear is almost 3.5 times lower than to the tension, biaxial component groups have from two up to four times smaller magnitude than the highest components and the lowest effective stiffness is reported for an interaction of shear and tension in different directions (Fig. 12). In spite of the coefficients corresponding to shear in the third model, all the effective tensor components linearly or nonlinearly decrease while increasing w with a relative reduction of 16% to almost 24% for w = 0.5. The interface defects are located in the interphase only, whose volume fraction in the composite is 6.5% and, therefore, such a reduction of the effective tensor is noticeable especially because it sharply increases with an increase in volume fraction of the interphase [13]. A difference in between the RVE computer models is the highest for the diagonal effective tensor components, especially for uniaxial tension, but the slope differs the most for shear deformations between the isotropic and anisotropic effective tensors.
The effective tensor components related to uniaxial tension and shear (Figs. 6, 7) generally decrease together with w in almost linear way and have the highest magnitude either for the third model (tension) or the first model (shear). These corresponding to the tension stay within 10% of each other, while the ones corresponding to the  Effective tensor components corresponding to biaxial tension for three spatial distributions of particles in the RVE shear are 3% more spread and much less convergent for the collected models. This is especially remarkable for the components of the third model, which behave in a totally different manner than the other ones quite contrary to the general trend. They start at 0.2 · 10 5 Pa, increase until a local maximum for w = 0.2 and then descend in a concave manner up to w = 0.5. Interestingly, these effective components cross the corresponding ones from the other models at w = 0.4 to become the highest above this volume fraction of interface defects.
The effective tensor components related to a biaxial tension (Fig. 9) and biaxial shear (Fig. 10) are decreasing with an increase of w in a manner close to the linear one and with total reduction from 20% to 26% for w = 0.5. These corresponding to biaxial tension are more than three times larger and start from a single magnitude of 3.42 · 10 6 Pa for w = 0 with a different slope for each model and almost perfectly coincide for each model. Interestingly, the results of the third model undergo a small fluctuation in the interval w ∈ 0.1..0.2 with an inflection point at approximately w=0.14. A situation is slightly different for the effective tensor components corresponding to biaxial shear-they start from two distinct values at w = 0, where C 45 is 22% higher than the other components and this difference decreases together with an increase in w; in this case, the various models give almost the same results for the entire range of interface defects parameter w.
The effective tensor components corresponding to two different stretches-tension and shear-in single direction (Fig. 11) and two different directions (Fig. 12) are also all decreasing with an increase in w and all have a very close result for all the coefficients coming from each model with a difference much smaller than a single percent. The results computed for extension in a single direction are twice as high as these from the  Fig. 11 Effective tensor components corresponding to an axial tension combined with shear in the same direction for three spatial distributions of particles in the RVE different directions and are 30% less affected by an increase in the parameter w; they also exhibit larger differences in between three given computational models. The effective tensor components obtained for uniform stretches in different direction converge at w = 0 to almost the same value of 1.71 · 10 6 Pa and comprise a componentC 26 , which for the second model diverges from the others for larger values of the parameter w.
Next, we compute the two engineering constants, i.e., Young's modulus E i and Poisson ratio ν i for all three models and plot them in Figs. 13 and 14, respectively. They are computed here mainly for the purpose of comparison and are not intended for representation of anisotropic response of the composite in full. These constants are both computed accordingly to an orthotropic solid in the first and the second model and the cubic solid for the third one. A set of relations for an orthotropic solid is the following one: It is composed of three various Young's moduli and six Poisson ratios. Please note that the tensional-shearing biaxial terms of the effective stiffness tensor (group 5 and 6 from Fig. 3) are disregarded in this computation, because for an orthotropic solid they are null. .
Plots in Figs. 13 and 14 additionally propose an average value for all the respective components E i and ν i to simplify the response into two elastic, well-known coefficients; they are all plotted in black color. Figure 13 firstly shows that Young moduli of the third model are considerably higher than for the anisotropic models, and it is less affected by an increase in w. It is also clearly visible that anisotropy in particle placement and clustering does not cause a high spread of this modulus in different directions. This is not true for Poisson ratios, visualized on Fig. 14, which show a considerable difference in magnitude for different directions, especially in the second model (with particle clustering). These ratios are generally higher for the anisotropic models, especially the 1st model, and decrease with w with a little higher rate to the third model. By this, the difference between the models lessens together with an increase in w. Next, we measure the anisotropy existent in the three composite models. This is done with use of three commonly used indexes, i.e., Zener anisotropy index A Z , Chung-Buessem anisotropy index A C and the universal elastic anisotropy index A U together with the new anisotropy measure, the tensor anisotropy index A T . Chung-Buessem and universal elastic anisotropy indexes are computed on the base of Reuss and Voigt bounds, calculated here in the following way: where K x and G x are the bulk and shear moduli for the composite according to Voigt and Reuss, K x¸G x denote the moduli of the three phases of the composite and f x means their volume fraction. The Chung-Buessem index is computed as [34] and the Universal elastic anisotropy index-with use of the below formula [35] The two next indexes are based solely on the basis of stiffness tensor (components) of the homogenized material. The first of these, Zener anisotropy index is defined in the following way [36]: The second one, called tensor anisotropy index A T , is derived directly from tensor coefficients and composed of following positively definite parts: It is introduced to leverage the limitations of A Z to the cubic tensors and defined as a sum of anisotropy coming from coefficients of group 1-3 A I (see Fig. 2) that are always present in a stiffness tensor and these coming from coefficients of group 4-6 A A -null for the isotropic, cubic, transversely isotropic and orthotropic tensors and positive otherwise. In this way, A I is enough detailed for measure of most tensors with symmetries and an additional term A A is meaningful solely for single symmetry (triclinic) and monoclinic tensors. A I is composed of the modified Zener index and a variance index A I = A I,z + A I,cov , where and where C eff Gi represents the ith group of tensor coefficients (see Fig. 2). A I,z defines a mean Zener index for all the three directions of the material, which simplifies to the regular Zener index A Z for isotropic and cubic tensor. A I,cov accounts for differences within the groups of coefficients 1-3, defining a sum of its coefficients of variation α. It is null for isotropic and cubic stiffness tensor and positive otherwise. The A A is defined as and it is intended for quantification of anisotropy coming from existence of coefficients that are null in the isotropic case together with their variation within these groups of coefficients (groups 4-6). The first part A A,m takes into consideration the expected value from these coefficients E C eff G4 , C eff G5 , C eff G6 that is multiplied by quotient n/12 relating number of nonzero coefficients n to its maximum number − 12. The second part A A,cov accounts for directional non-uniformity that causes different response for the same action in three axes and it is defined as a sum of variances for groups 4-6. Both parts of A A are divided by the expected value of coefficients from the first three groups E C eff G1 , C eff G2 , C eff G3 so that the additional anisotropy caused by these terms is referenced to the principal coefficients (principal response of the continuum). This is done to reflect the case where the additional effects from groups 4-6 are relatively small in comparison to the main response, while quotient n 12 accounts for the case when only several coefficients are not null (for example in the case of triclinic tensors) not to overestimate their effect. A full formula describing this index has been introduced as Such a coefficient could be computed for all types of elastic stiffness tensors and, unlike all remaining indexes, it takes into consideration all the stiffness tensor coefficients that are not null in the specific tensor. The main advantages of this measure can be listed as follows: single valued, modular-complexity of calculations increases together with an increase of complexity of the analyzed stiffness tensor, derived directly from tensor properties, universal-applicable for all kinds of tensors (not limited to cubic, orthotropic or other types of tensor), doesn't require additional calculations of the Reuss-Voigt bounds, versatile-applicable not only for anisotropy in the stiffness tensor but also to other tensors, can be applied for homogenized or not homogenized materials alike and, unlike any other measure, it takes into consideration all components of this tensor that are not null, allows quantification of anisotropy of arbitrary genesis, for example those coming from spatial composition of the reinforcement or non-uniform voids, for which the other measures are not perfect (see A U and A C on Fig. 15, whose measure is exactly the same regardless the spatial composition of reinforcement). The minimum value of this measure is 1.0 and it is equivalent to a full isotropy of the tensor. The above measures of anisotropy are used for the three models of the composite and presented in Fig. 15 in relation to the volume fraction of voids w. These include the Zener anisotropy index A Z , Chung-Buessem anisotropy index A C , Universal anisotropy index A U and the new tensorial anisotropy index A T that are differentiated by symbols on the graph. The A Z and A T are computed for all three models marked by colors, while the A C and A U presented only once as they return the same results for all models. This equality of results comes directly from the formulas of A C and A U that are based on the Reuss-Voigt bounds, independent of the spatial composition of reinforcement in the RVE. We also note the high magnitude of A U that comes from a considerable difference of bulk and shear moduli of the Reuss-Voigt bounds, whose quotients are calculated in this measure. For better readability of the results on Fig. 15 we present a quotient A U /1000, not A U . Looking on this figure, we first report that the A C and A U are practically independent of w, while the remaining two methods, based directly on the effective tensor, report a decrease in anisotropy together with an increase of w (they become closer to isotropic). This finding is quite important not only for the linear regime of the composite, it is also interesting from the view of fracture mechanics. This is because micro-damage in such composites initiates from the interphase and, quite similarly to the voids, it weakens connection between the particle and the matrix effectively making such anisotropic composite closer to isotropic. Magnitudes of these four measures are wide spread; A U ∼ = 1490 is the highest one, A C ∼ = 0.99-the lowest and A Z , A T ∈ 1.0, 1.8 . These two last measures (A Z and A T ) also report major differences in anisotropy for the three models-the 1st and the 2nd both have a considerable anisotropy, while the 3rd one is perfectly isotropic-A Z = A T = 1.0 with only a small oscillation for several data points (presumably numerical errors). This result follows the spatial placement of the reinforcement, isotropic for the 3rd model and anisotropic for the two remaining ones (see Figs. 2,8).The tensorial anisotropy index reports a higher anisotropy for the first two models than the A Z , because it takes into consideration also the additional terms for complex tensors. This difference is quite high because the considered tensor is monoclinic, has axial variation of properties and all 21 tensor coefficients present. The difference of anisotropy in the first two models is moderate, marginally dependent on w and recognized in similar manner in both, A Z and A T ; Last but not least, Fig. 15 shows that clustering decreases a little the degree of anisotropy (magnitude smaller for the 1st model than for the 2nd one) but it is also less affected by an increase of w-it decreases a little slower than for the model without clustering, which is evident for, A Z and A T .
Further, sensitivity coefficients S C eff αβ of all effective tensor components (Figs. 16,17,18,19,20,21) are calculated from the polynomial-based response functions with use of a direct algebraic differentiation in MAPLE. They are defined in the neighborhood of one of the expected values of the volume fraction of interface defects E (w) ∈ [0.1, 0.2, ..., 0.5] and are non-dimensional. Therefore, they represent a set of local sensitivities of the respective effective tensor components. This is the reason why, in some specific cases, they do not catch perfectly the general tendency of this components. This takes place for sensitivities corresponding to the pure shear of the 3rd model and indicates some local fluctuation of majority of the others around the expected value of the parameter w. Sensitivity coefficients comprise 55 independent derivations of discrete data points and are composed from 5 response functions for each legend item. These, in turn, serve for conversion of the discrete results of the effective tensor components C eff αβ into their continuous representation in relation with the  It is quite expected after previous considerations that almost all the sensitivities are negative for an entire considered spectrum of the volume fraction of voids w, decrease together with an increase of w and start close to but below 0 (typically at S C eff αβ = −0.05). They decrease until 0.2 ÷ 0.35 for w = 0.5, are very close for the anisotropic models and always smaller for the isotropic composite model. This means that all the Fig. 17 Sensitivity of effective tensor components corresponding to a pure shear of the RVE w.r.t. the interface defects ratio w effective tensor components increase together with an increasing w and preserve a very close slope (magnitude of sensitivity is very low). This is not entirely true for the third model, which returns some local effects for several effective components. All the results are smooth and continuous with an exception to the ones of the 3rd (isotropic) model that are either minimally (as for Fig. 16) or heavily (Fig. 17) fluctuating together with increasing parameter w. Interestingly, this high fluctuation seems to concentrate around a certain magnitude of S C eff αβ = 0.05 and is well reflected by the dispersion of results of the effective tensor component in Fig. 7. This all means that for the anisotropy the effective tensor components are more affected by the interface defects and this relation with w has less local effects than for simple isotropic model. Sensitivities corresponding to uniaxial tension and shear (Figs. 16,17) firstly show almost a perfect agreement of the results coming from the first and the second model and dispersion of these coming from the third model for each C (eff) . This fluctuation indicates a local character of the sensitivities coming from the third model and since that they require a special computational treatment. The remaining two models also have some small fluctuation for each E (w) coming from the same source, but in all cases it does not diminish clarity of the results. A magnitude of sensitivity for these models is up to 25% smaller for the effective tensor components corresponding to shear. A situation is entirely different for the 3rd model not only due to the positive values, but also because of the fact that the sensitivity magnitude corresponding to uniaxial tension is up to several times smaller. This positive value is reflected by an initial increase of the effective tensor component C eff 44 -the sole positive response of the effective tensor to an increase in the volume of interface defects.
Sensitivities corresponding to biaxial tension and shear (Figs. 18, 19) are a little less convergent for the anisotropic models, and the third model returns very similar sensitivity to the remaining ones; these coefficients have a very similar magnitude, the same sign and a close slope. They reach 0.34 for w = 0.54 and are almost linear with respect to the parameter w. These response functions have minimally more local character than for their uniaxial counterparts, especially for the first two models, while the third model has almost no fluctuations except for E (w) = 0.1, which returns a considerable higher value than the others. Sensitivity coefficient S C eff 46 is several percent smaller than the remaining coefficients corresponding to shear for both available models-this difference is a little larger for the second model and increases together with an increase in w. This is in a perfect accordance with Figs. 10 and 11, where the slopes of all the effective tensor components are very close. Even these small differences in sensitivity do not have a high impact on their behavior in relation to the interface defects ratio w; such a difference is not remarkable for the sensitivities corresponding to a biaxial tension, which differ with about 1% of each other.
Sensitivities coming from a simultaneous action of tension and shear in the same or different direction (Figs. 20, 21) are all similarly affected by the volume fraction of the interface defects. They all decrease, have almost uniform magnitude and a negative sign for the respective models and effective tensor coefficients. The former ones (Fig. 20) start approximately at S C eff αβ = − 0.03 and decrease up to S C eff αβ = − 0.25, while the latter ones (Fig. 21)-at S C eff αβ = −0.05 and S C eff αβ = − 0.31. The only divergence occurs for the 1st model and E (w) = 0.2, where some coefficients locally exhibit a sharp increase or decrease together with w. The 1st model returns a little bit larger results, but this difference is much below a single percent and thus negligible. All the sensitivity coefficients computed here feature a small curvature coming from their Fig. 20 Sensitivity of effective tensor components corresponding to an axial tension combined with shear in the same direction for the RVE w.r.t. the interface defects ratio w Fig. 21 Sensitivity of effective tensor components corresponding to an axial tension combined with shear in different directions for the RVE w.r.t. the interface defects ratio w local character, but this does not prevent us to discern a clear downward and almost linear relation of these coefficients with respect to the parameter w.

Probabilistic computations using the ISFEM
The most important part of this analysis is determination of statistical disorder of the anisotropic and isotropic effective tensor to random variability in the volume fraction of interface defects w, which is done w.r.t. the coefficient of variation of this random parameter α (w) and for its increasing expected value E (w) in the range of w ∈ [0.10, 0.20, . . . , 0.50]. The series of FEM experiments have been carried out for this purpose using the WLSM technique to determine the polynomial response function and statistically optimize its order. The resulting polynomials are integrated together with the Gaussian PDF in a semi-analytical method, inserted into the 2nd-and 4th-order Taylor expansions in stochastic perturbation-based approach and finally subjected to a discrete Monte Carlo sampling with 250.000 trials for an increasing α (w) . The response function serves here as a conversion of the discrete results of the FEM computations for a continuous representation of the effective tensor components as the functions of the interface defects volume fraction (in the interphase). An optimization process of this polynomial function consists of a simultaneous maximization of the correlation and minimization of both the variance and RMS error of the WLSM. This ensures the best fitting of the polynomial to the discrete computation results and such an optimization has been made separately for each tensor component in all models and each E (w) for full polynomials of orders from the first until the tenth. Optimization and WLSM computations have been programmed in algebra system MAPLE for the Dirac-type weight distribution w e ∈ [1, 1, 1, 1, 1, 10, 1, 1, 1, 1, 1], which gives the major weight for the expected value being equal to a sum of the weights associated to the rest of discrete values in its 10% neighborhood. It is necessary to underline that the computational effort in semi-analytical considerations and perturbation-based approaches (2nd and 4th order) are negligible when compared with this corresponding to the Monte Carlo simulation (MCS).
Probabilistic considerations are here twofold. First of all, we compare the influence of the input coefficient of variation of interface defects fraction α (w) on the resulting probabilistic characteristics of the effective tensor components. These are plotted consecutively on the left side of Figs. 22-45, where the colors define the expected value of volume fraction of interface defects E (w), which the computations are made for and are marked in the legend as w = 0.x. Secondly, we examine coincidence of the probabilistic characteristics made by the three most common probabilistic methods -the semi-analytical method (AM), the Stochastic Perturbation method of the 2nd order (SPT 2nd), as well as 4th order (SPT 4th) and the Monte Carlo simulation (MCS). These are placed on the right side of all these graphs and presented for one expected value of interface defects fraction E(w) = 0.3 only, which is done for brevity of the presentation. Final computations include the first four probabilistic characteristics of a chosen effective component C (eff) from each group and these are the expected values (Figs. 22, 23, 24, 25, 26, 27), coefficients of variations (Figs. 28, 29, 30, 31, 32, 33), skewness  (Figs. 40, 41, 42, 43, 44, 45). All these coefficients and moments are non-dimensional, with an exception to the expected value which, similarly to the effective stiffness tensor, has a unit of Pascal. All the trends obtained with use of the ISFEM demonstrate smooth sufficiently variations with some small or a complete lack of local fluctuation at all with an exception to the skewness and kurtoses, which in some cases have a local vertical asymptotes. However a coincidence of the expectations is almost perfect and it concerns the 4th -order method after 4th-order response polynomial-the SPT of the 2nd order diverges for most probabilistic characteristics where α (w) > 0.05. The coefficients of variations are not so perfect, but still the AM and MCS are in almost perfect agreement, the SPT being a little less precise for the input coefficients of variation when α (w) > 0.08. The SPT of the 4th order gives definitely a better result than the one of the 2nd order and a higher order would be necessary for a better match with the other probabilistic methods. The results are a little dispersed for the higher probabilistic characteristics of effective tensor, where not only the SPT but also the remaining two methods lose their perfect correlation. The expected value of some effective tensor components has a little bit smaller (about 3%) dependence on the random volume fraction of interface defects than their deterministic counterparts. Further, we observe that the coefficients of variation α C (eff) αβ in almost all the cases are 2-3 times smaller than the input coefficient of variation α (w) . This, however, does not mean that the increase in Gaussian dispersion of w is irrelevant for the magnitude of E C eff αβ , because it not only changes this expected value up to 15% for α = 0.15, but also result in a considerable dispersion of the random effective tensor. This dispersion cannot be bounded by the 3σ rule, because the output is not Gaussian after initial theoretical considerations at least.
Expected values of the effective tensor components (Figs. 22,23,24,25,26,27) are decreasing all together with an increase in E (w), are moderately affected by the α (w) and show a considerable differences in between the three RVE models. It is seen not only in the magnitude of the expected values of E C eff αβ , but also in the slope and curvature of relation with respect to α (w) , which can be either decreasing and concave as well as increasing and convex. The largest expected values result either from the third model or from the first model; the first model return also extreme values of effective tensor. The differences, however, never exceed 11% and for all but one representative component they are in range of 1-4%. Maximum reduction of E C eff αβ for E (w) = 0.5 goes up to 20% for E C eff 56 and anisotropic models, while the minimum one is 4% and occurs for E C eff 55 only for the third model. This is a little smaller than the highest decrease of the C eff αβ and implies a small reduction of the dependence of the effective tensor on the volume fraction of interface defects when considering its Gaussian dispersion. A typical reduction of E C eff αβ for E (w) = 0.1 is approximately 3% and is more than two times smaller for the isotropic model than for the anisotropic cases. On the other hand, isotropic model is typically several times more affected by α (w) and is the only one, where the expected values calculated for different (and especially adjacent) E (w) cross each other for α (w) > 0.1. One more interesting observation is that expectation E C eff 55 in the 3rd model not necessarily decreases together with an increase in E (w) and this phenomenon has been reported before in Fig. 7. Similarly to the effective tensor components, the most transparent differences in between probabilistic methods are noticed for the uniaxial counterparts of the expected value E C eff αβ . A comparison of different probabilistic methods show their perfect agreement for all α (w) with an exception to the SPT of the 2nd order, which starts to diverge at α (w) > 0.06 for all the components. This clearly indicates that the 2nd order is rather ineffective and the 4th-order expansion is needed. Expected values of the isotropic composite are much closer to each other-they are much less suscep- tible to w than for the anisotropic models. Moreover, their values are either higher than (E C eff 22 , Fig. 22 and E C eff 12 , Fig. 24) or a little beneath the middle of mean (E C eff 55 , Fig. 23) of the anisotropic expectations. For other groups of coefficients isotropic counterparts are unavailable.
The coefficients of variation of the effective tensor components (Figs. 28,29,30,31,32,33) are all positive and increasing in a convex manner together with an increasing α (w) . They are always lower than the input coefficient of variation and increase their magnitude together with an increase of E (w) . The only exception is noticed for C eff 55 corresponding to the 3rd model, which reaches α C eff αβ = 0.26 for E (w) = 0.5. The results obtained for this model are either much higher than α C eff 55 , comparable to α C eff 22 or a little smaller than α C eff 12 these coming from the anisotropic models. The anisotropic models generally return the results very close to each other, but their coincidence decreases together with both, an increasing α (w) and E (w). There neither exists a single model with the highest output coefficient of variation, nor with the lowest, but the third model has the biggest scattering of the results for the given E (w) reaching both the highest and the lowest values. All the coefficients of variation should and begin from 0 and increase up to 0.005 for E (w) = 0.1 and 0.05 for E (w) = 0.5 at α (w) = 0.15. This difference is with about ten times, but the largest output coefficients of variation are almost five times smaller than the input ones. This means that the output effective tensor coefficients are generally less distant from each other than the input volume fraction of interface defects. This difference, however, increases together with an increasing volume fraction of the interface defects. The results of three probabilistic methods are close to each other, but the SPT no longer keeps a perfect match for α (w) > 0.08 (Figs. 28b, 29b, 30b, 31b, 32b, 33b)-especially the SPT of 2nd order. Even the results of the 4thorder SPT are a little bit lower than these obtained from the MCS and AM, which still are highly coincident. This enables to recommend further usage of the SPT of 4th order in analysis of the coefficient of variation relative to E (w) with a very small error and, consecutively, such a method is chosen to prepare the results shown in the left parts of Figs. 28, 29, 30, 31, 32 and 33. The isotropic model has generally higher coefficient of variation of the stiffness tensor than the anisotropic ones, with an exception to α C eff 22 , which stays a little lower. Skewness of the effective tensor components (Figs. 34,35,36,37,38,39) is remarkably influenced by the input coefficient of variation α (w). Its relation with this coefficient has, nonetheless, complex character as for the lower probabilistic moments-it either has a local maximum at α (w) = 0.25, as for β C eff 55 , or increases in a convex or concave manner. It is neither predominantly positive nor negative for each effective tensor coefficient but, in most cases, it has the same sign for each model. The skewness computed for the third model is usually the largest, while the one corresponding to the first and second models has always very similar magnitude (although not the same sign). There exists no strict relation of the skewness with E (w) but it is close to 0 for the anisotropic models and α (w) < 0.7; only then it starts to sharply increase or decrease and reaches values of β C eff αβ between − 6 and 9. It needs to be underlined that a coincidence of the three chosen probabilistic methods for skewness has rather conditional-they all return the same result for the input coefficient of variation up to α (w) = 0.04, but diverge for higher α (w) . This is especially true for the SPT, which always underestimates the skewness, but, with an exception to β C eff 12 , it returns the same sign as the other methods. It is quite important, however, that also the two remaining probabilistic methods do not coincide well-the MCS returns larger skewness for almost all the effective tensor coefficients. This is true for all the RVE models, but the first one is here the one with the best agreement of all computational methods and the third one-with the worst. Therefore, unlike for the other moments, the semi-analytical method (AM) is chosen for representation of the influence of E (w) on all the effective tensor components C (eff) . Skewness of the isotropic model has a higher magnitude than the ones computed for the anisotropic ones. It is also a little more convex in relation to α (w) and, in contrast to the skewness coming from anisotropic models, it has an apparent local maximum for β C eff 55 obtained for α (w) = 0.02 and exists for all mean values of volume fractions of the defects. Kurtosis of the effective tensor components (Figs. 40, 41, 42, 43, 44, 45) are all positive, start from zero and usually reach a considerable magnitude for high α (w)-even higher than the skewness. They are either convex and always increasing or concave and with an apparent limit of 3; this feature is neither model nor component specific. The kurtosis typically increases together with an increase of E (w) , but this is not always true-sometimes it increases for small E (w) and then decreases from E (w) = 0.3 or E (w) = 0.4. It is not clear, which of the models return the highest or the lowest results, but the isotropic model usually has a higher magnitude then the anisotropic ones for low α (w) and lower for α (w) > 0.11. Similarly to the skewness, results of the three probabilistic methods are not in a perfect agreement for the kurtosis. Despite the use of 250.000 trials, the MCS is not more precise than the SPT for the kurtosis and they both are sometimes distant from the AM. On the other hand, this last method shows some discontinuities especially for α (w) close to 0, where stochastic problem becomes deterministic and which is quite typical and expected for the AM with a response function having polynomial relation with an input argument. The SPT is preferable for use in analysis of κ C eff αβ respective to the mean value of volume fraction of interface defects E (w) principally, because it is the quickest one and additionally because it gives a continuous representation of the kurtosis without any discontinuities. The isotropic model has either a similar (κ C eff 22 ) or higher kurtosis (κ C eff 55 and κ C eff 12 ) than the anisotropic models. It also tend to have some apparent maxima of kurtosis especially for κ C eff 55 , which are not reported for anisotropic models at all.
The random effective tensor for all the considered models cannot have Gaussian distribution and this observation is consistent with the previous analytical derivations. This is because both skewness and kurtosis remarkably differ from 0 and typically increases together with an input uncertainty. Its variation, however, is much smaller than the input for almost all the effective tensor components. An uncertainty in w affects the isotropic effective tensor a little more than the anisotropic one. This is remarkable especially for the expected values, which are much more affected by the coefficient α (w) than these resulting from the anisotropic models and in the coefficients of variation, also usually higher for the isotropic case. These differences are increasing together with an additional increase in expected value of the interface defects fraction E (w) and its uncertainty in the composite.

Conclusions
Probabilistic homogenization of the particle-reinforced composite with stochastic interface defects and anisotropic distribution of the particles has been demonstrated in this work by using of three independent stochastic computer techniques. The original semi-spherical defects with randomly defined radius having Gaussian PDF have been replaced by the isotropic interphase of constant thickness having randomized material characteristics. One of the most important conclusion of homogenization of such a composite is that effective elasticity tensor components cannot be Gaussian in this case and this is not for sure an effect of the particles distribution anisotropy. It is important that this fact has been independently proven analytically and also verified numerically. This paper also proposes a new measure of anisotropy called the tensorial anisotropy index, which is able to capture the difference in anisotropy caused by spatial placement of its phases and answers the need for anisotropy index effective beyond the cubic tensors (Fig. 15).
There is no doubt that overall computational cost of anisotropy effect modeling is much larger than for the isotropic effective elasticity tensor because instead of two really independent components, one must calculate 21 components, which is especially challenging during stochastic simulations according to the unweighted statistical schemes. FEM mesh and its generation require at least the few particles to be interconnected by non-uniform meshes parts for the matrix region and it requires more than 10 5 additional hexahedral or tetrahedral elements into the entire RVE and separate meshing around any of 27 particles; it would be recommended further to treat a single carbon particle as the super-element in ABAQUS model of the multi-particle systems static or dynamic equilibrium as they are nearly non-deformable due to the large contrast of elastic modulus in addition to the matrix. This effort is compensated by a more precise determination of the stiffness tensor which, for anisotropic placement of particles, also proves anisotropic and ensures a much better prediction of the internal stress intensity vital in fracture mechanics and durability prediction. Let us note that in our case the Zener coefficient of anisotropic model reaches 1.12, while the tensorial anisotropy index goes up to 1.78 and they are dependent on the volume fraction of voids, whose increase decreases stiffness anisotropy. Let us note that the measures of anisotropy based on the Voigt-Reuss bounds are ineffective in determination of an influence of spatial placement of particles on the anisotropy in the effective (homogenized) stiffness tensor, as they return exactly the same result for all the three models.
Differences between the stiffness tensor of the isotropic and anisotropic RVEs are determined not only for the deterministic case, it is also evident for its probabilistic moments and coefficients. The expectations of isotropic composite are much less dependent on the expected value of volume fraction and are generally higher for the isotropic RVE, coefficients of variation are also a little higher, skewness has apparent local maxima and kurtosis-global ones, both not existent for the anisotropic models. When comparing the engineering constants of these models, i.e., Young's modulus and Poisson ratio anisotropy decreases the first one but increases the second. A weakening of the Young's modulus in anisotropic models increases together with an increase in volume fraction of voids in the interphase and so does the strengthening effect of the Poisson ratio. It is also worth to mention that an effect of anisotropy on the effective stiffness tensor is rather global, which means that it affects entire groups of coefficients and does not provide a very significant difference between the stiffness in different directions, at least in the linear regime. Please note that reduction of 21 coefficients to 6-each coming from different group of coefficients-is still a lossy process.
Further one needs to notice that the optimal response polynomials of C (eff) αβ in addition to the interface defects ratio w for both isotropic and anisotropic RVEs homogenization are the family of the same 4th-order polynomials (according to the statistically optimized WLSM technique). It was the basis to derive fundamental equations for probabilistic moments of the effective elasticity tensor in addition to the basic characteristics of uncertain parameter w. The effective tensor in anisotropic case generally is more sensitive to this parameter fluctuations than in the isotropic model. The expectations E C (eff) 11 and E C (eff) 12 appear to be always larger for isotropic model than in the anisotropic one, while E C (eff) 44 shows an inverse tendency, at least up to the limit value w=0.40. Homogenization of the given composite remarkably damps initial uncertainty in interface defects and it follows a local character of this variable within the RVE (even after probabilistic averaging throughout the interphase), so that probabilistic entropy as a universal measure of statistical disorder of the RVE decreases with the few times during homogenization process. As it was expected after previous theoretical [14,15] and numerical studies [7,8] a presence of the interface defects and their uncertainty significantly reduces the mean values of the effective elasticity tensor also for the particle-reinforced composites.
Probabilistic homogenization method displayed in this work may also benefit with determination of the cross-correlations of the effective elasticity tensor components thanks to the semi-analytically determined polynomial responses of these components to Gaussian w. The 3D FEM model of the particle-reinforced composite with anisotropic distribution of the particles prepared for this study may further serve with almost no modifications to stochastic simulation of uni-and bidirectional stretches of the RVE, whose first trial has already been performed [37]. Assuming that the variable R exhibits a Gaussian distribution with the first two moments equal to E[R] and Var(R), one can derive by an integration from the definitions of these two moments or, alternatively, using the complex characteristic function of this distribution, as It is easy to prove that the induced random variable w cannot have Gaussian distribution in this case, because lim w→0 w − 2 3 = +∞, so that the probability density of w itself is unbounded. Recalling further the Gaussian upper bound on the interphase thickness calculated from the interface defects radius as one can derive the first two probabilistic moments of the variable w as Var (w) = 3 n 2 Var 2 (R) 3E 4 [R] + 12E 2 [R] Var(R) + 5Var 2 (R) and substitute them into further computations devoted to homogenization at the microscale. One needs to recall the fact that the interface defects radius is really distributed according to the truncated Gaussian distribution, where at least negative values are prohibited as physically disabled. Upper and lower limits for probabilistic integrals describing expectation and variance of the parameter w are derived analytically using the same three-sigma rule as for the interphase thickness, i.e., This leads to the following algebraic formula: where a denominator D equals to It should be mentioned that rational parts of the coefficients included in the above equations (29)(30) correspond to the analogous coefficients in the expanded version of equation (26); these coefficients of course have no physical meaning. Moreover, the additional formula for the variance including truncation error can be derived in a similar way but has many times larger expression and so that it is postponed here. Therefore, we assess truncation modeling errors of the expectation and the variance of w by calculating relative errors for both functions. The following formulas do apply: Finally, we plot these errors in Figs. 46 and 47 for unit particle radius r =1 and ten interface defects n=10. As one may expect, both errors are positive despite the input parameters choice and ranges from the few promiles for the expectations up to the few percents − for the variances and obviously both equal 0 for deterministic case when α(R) vanishes. Truncation error in expectations increases together with this parameter and is totally insensitive to expectation of the interface defects mean radius. The variance error demonstrates the same insensitivity as for E[R], whereas sensitivity to α(R) is very huge close to the deterministic origin (throughout the interval α(R) ∈ [0.00, 0.01]) and then splashes and becomes almost linear. Positive values of both errors mean that the models without truncation effects exhibit larger expectations and variances of the parameter w, so that it is impossible to verify whether they under or overestimate the interphase effect on the RVE response in its stochastic finite element method analysis; nevertheless, it should not be simply postponed, especially for the non-Gaussian random variables.