Energy Fluctuations in the Homogenized Hyper-Elastic Particulate Composites with Stochastic Interface Defects

: The principle aim of this study is to analyze deformation energy of hyper-elastic particulate composites, which is the basis for their further probabilistic homogenization. These composites have some uncertain interface defects, which are modelled as small semi-spheres with random radius and with bases positioned on the particle-matrix interface. These defects are smeared into thin layer of the interphase surrounding the reinforcing particle introduced as the third component of this composite. Matrix properties are determined from the experimental tests of Laripur LPR 5020 High Density Polyurethane (HDPU). It is strengthened with the Carbon Black particles of spherical shape. The Arruda–Boyce potential has been selected for numerical experiments as ﬁtting the best stress-strain curves for the matrix behavior. A homogenization procedure is numerically implemented using the cubic Representative Volume Element (RVE). Spherical particle is located centrally, and computations of deformation energy probabilistic characteristics are carried out using the Iterative Stochastic Finite Element Method (ISFEM). This ISFEM is implemented in the algebra system MAPLE 2019 as dual approach based upon the stochastic perturbation method and, independently, upon a classical Monte-Carlo simulation, and uniform uniaxial deformations of this RVE are determined in the system ABAQUS and its 20-noded solid hexahedral ﬁnite elements. Computational experiments include initial deterministic numerical error analysis and the basic probabilistic characteristics, i.e., expectations, deviations, skewness and kurtosis of the deformation energy. They are performed for various expected values of the defects volume fraction. We analyze numerically (1) if randomness of homogenized deformation energy can correspond to the normal distribution, (2) how variability of the interface defects volume fraction a ﬀ ects the deterministic and stochastic characteristics of composite deformation energy and (3) whether the stochastic perturbation method is e ﬃ cient in deformation energy computations (and in FEM analysis) of hyper-elastic media.


Introduction
Interface defects are the very important engineering problem in composite materials. This is because they are the result of many manufacturing processes, they remarkably affect their mechanical (and coupled) response and occupy relatively small total volume fraction of the entire composite. They affect failure [1,2], reliability [3], durability [4] and thermal conductivity [5] of the composite and can be included in numerical analysis as geometrical imperfections [6]. Alternatively, interface defects can be represented as imperfect interfaces, which commonly serve for a realistic prediction of composites behavior [7,8]. They can also be adjusted to fine tune the characteristics of composites for specific needs [1,9]. Defective composites are investigated with an additional thin layer between its initial constituents [8,[10][11][12][13] defined as the interphase. This model follows some manufacturing processes relevant to the polymeric composites where glassy phases may occur around the reinforcing fibers of particles due to thermal processing of such a composite. A very similar interphase can be noticed in rubber-based mixtures, where the occluded rubber forms some kind of the thin interphase about the reinforcing carbon particles. Other numerical approaches to the interface defects include introduction of a system of springs [14] or addition of specific contact finite elements [15]. As it has been documented, this interphase is the weakest link in composite structures and significantly affects effective (homogenized) material properties of multiple composites irrespective to their degree of anisotropy and also outside of the deterministic context. They affect probabilistic and stochastic characteristics of these composites [16,17]. An increase of probabilistic dispersion of defects volume fraction results in a considerable uncertainty in mechanical properties of entire composite and an increase of their amount typically causes degradation of these properties. An influence of the interphase could be either negative, in the case when the bond between the main constituents is weak or defective [18], or positive in the presence of chemical bonds between the principal phases [19,20], for example in bound rubber. Early approach to homogenization of composites included a pure analytical apparatus. It was strictly limited in application to linear elastic, reversible and isotropic constituents and also required a simple pattern of periodicity of the Representative Volume Element (RVE) [21][22][23]. Contemporary studies use a multi-scale approach rather to represent a nonlinear response for composites, elasto-plastic, visco-elasto-plastic and also the hyper-elastic one [24][25][26][27]. In such studies, deformation energy is necessary for determination of nonlinear model parameters but only few of them [27,28] focus on it directly. Nevertheless, they do not address the fundamental question of how the interface defects affect deformation process uncertainty of the composite and are not based on experimental tests for a specific material, which is done in this study.
The principle aim of this work is to determine how the initial stochastic interface defects affect the deformation process of a composite under classical uniaxial tension experiment; a dual probabilistic approach to the hyper-elastic problem is proposed here to achieve this goal. This is done to validate the Iterative Stochastic Finite Element Method (ISFEM)computation against the classical Monte-Carlo simulation [29] and provide an efficient alternative offering a reduced computational time and requiring less computer resources. The ISFEM is based on the iterative stochastic perturbation approach and is augmented with a supplementary FEM implementation. The deformation energy is determined with use of the homogenization technique and computed with the Finite Element Method (FEM) based on displacements. A unitary cubic Representative Volume Element (RVE) is used. The matrix of this composite is applied after laboratory tests of Laripur LPR 5020 and the particles are modeled as Carbon Black Fullerenes F60.
Stochastic perturbation approach is based on a cyclic FEM usage with a fluctuating parameter for numerical recovery of the polynomial response function relating the resulting deformation energy with the given uncertain initial parameter, i.e., interface defects volume fraction. The input parameter is included in a polynomial form directly in the effective strain energy of the composite that is based on Arruda-Boyce potential or, alternatively, this relation is represented as a bivariate polynomial. It is used for calculation of the first four probabilistic characteristics of the structural response using stochastic perturbation approach based upon the general order Taylor expansion.

Methods
Let us introduce a continuous and heterogeneous solid body in 3D Euclidean space Ω visualized on Figure 1. It describes a particle-reinforced composite and contains three distinct phases-matrix Ω m , interphase Ω i and particles Ω p and so that Energies 2020, 13, 2011

of 16
The matrix and an interphase are both hyper-elastic and particles linear elastic. It is further assumed that a contact in-between matrix and interphase as well as in-between the particle and the interphase too is perfect and also load-independent. Further, it is mandatory to assume that interphase and the particle does not intersect the outer surfaces of the RVE ∂Ω , so that ∅ = Ω ∂ ∩ Ω p . Total defects volume in the given RVE can be expressed via summation of all the defects contributions as follows is approximated with semi-spherical shape penetrating the matrix and a geometrical center of such a defect is each time located on the particle surface. Each defect radius R is assumed to be statistically dispersed according to the normal distribution and is defined by its expected value E[R] and variance Var(R). Next, such a defected interface is replaced with the interphase (as artificial layer around the particle) having constant thickness around the particle, which includes all the interface defects and small portions of the matrix surrounding them. This thickness also serves as the upper limit of the defects radii that is limited by a three sigma rule where ∆ is a thickness of hyper-elastic interphase which surrounds a particle. All material characteristics of the interphase necessary for further numerical simulation are simply calculated as volumetric average of the matrix and the defects contributions. Similarly, the effective stress of the interphase is recovered using the rule of mixtures connected to volumetric contributions of the matrix and the defects contained in Ω where denotes the stress intensity in the matrix contained in the interphase and stress intensity in the defects. The first two probabilistic moments of σ int are further represented as The matrix and an interphase are both hyper-elastic and particles linear elastic. It is further assumed that a contact in-between matrix and interphase as well as in-between the particle and the interphase too is perfect and also load-independent. Further, it is mandatory to assume that interphase and the particle does not intersect the outer surfaces of the RVE ∂Ω, so that Ω p ∩ ∂Ω = ∅.
Total defects volume in the given RVE can be expressed via summation of all the defects contributions as follows A single interface defect Ω d(i) is approximated with semi-spherical shape penetrating the matrix and a geometrical center of such a defect is each time located on the particle surface. Each defect radius R is assumed to be statistically dispersed according to the normal distribution and is defined by its expected value E[R] and variance Var(R). Next, such a defected interface is replaced with the interphase (as artificial layer around the particle) having constant thickness around the particle, which includes all the interface defects and small portions of the matrix surrounding them. This thickness also serves as the upper limit of the defects radii that is limited by a three sigma rule where ∆ is a thickness of hyper-elastic interphase which surrounds a particle. All material characteristics of the interphase necessary for further numerical simulation are simply calculated as volumetric average of the matrix and the defects contributions. Similarly, the effective stress of the interphase σ int is recovered using the rule of mixtures connected to volumetric contributions of the matrix and the defects contained in Ω i where σ m denotes the stress intensity in the matrix contained in the interphase and σ d stress intensity in the defects. The first two probabilistic moments of σ int are further represented as Let us next assume that defects do not provide a contribution for the stress σ d (ε i ) = 0 and introduce a parameter w denoting their volume fraction, so that (6) and let us assume that this parameter is independent of the strain level. Its characteristics, i.e., expectation E[w] and variance Var[w] are determined analytically from E[R] and Var(R) using the basic definitions of probability theory; some details may be found in Appendix of the reference [16]. Multi-scale homogenization scheme is proposed to model their impact on the composite in its macro-scale and it consists of three steps. Firstly, the interphase is replaced with the new artificial material surrounding all the particles. Secondly, the interface defects are probabilistically averaged throughout the corresponding interphases. Finally, the composite is probabilistically homogenized together with all these interphases. A micro-geometry of this composite has been schematically given in Figure 2. Deformation gradient of hyper-elastic matrix is defined as F = ∂x/∂X, where x represents the spatial configuration and X describes the reference configuration and deformation energy accumulated in the matrix can be described as where F T F is the right Cauchy-Green deformation tensor. Let us next assume that defects do not provide a contribution for the stress ( ) = 0 and introduce a parameter denoting their volume fraction, so that = + = (6) and let us assume that this parameter is independent of the strain level. Its characteristics, i.e., expectation [ ] and variance [ ] are determined analytically from E[R] and Var(R) using the basic definitions of probability theory; some details may be found in Appendix of the reference [16]. Multi-scale homogenization scheme is proposed to model their impact on the composite in its macroscale and it consists of three steps. Firstly, the interphase is replaced with the new artificial material surrounding all the particles. Secondly, the interface defects are probabilistically averaged throughout the corresponding interphases. Finally, the composite is probabilistically homogenized together with all these interphases. A micro-geometry of this composite has been schematically given in Figure 2. Deformation gradient of hyper-elastic matrix is defined as = / , where represents the spatial configuration and describes the reference configuration and deformation energy accumulated in the matrix can be described as where is the right Cauchy-Green deformation tensor. Let us assume that this matrix is incompressible, its behavior depends upon the first two invariants only and that it can be efficiently represented by the Arruda-Boyce model [30]. Then, its energy can be calculated as where 1 is an additional material parameter, ̂= 1 2 , stands here for the locking stretch of the polymer chain networks, while define the additional coefficients (given in Table 1). The effective strain energy of homogenized composite depends on the invariants and interface defects volume fraction. It could further be divided into three independent contributions-of the matrix, of the interphase and of the particle, independently. There holds Let us assume that this matrix is incompressible, its behavior depends upon the first two invariants only and that it can be efficiently represented by the Arruda-Boyce model [30]. Then, its energy can be calculated as where C m 1 is an additional material parameter,β = 1 , λ m stands here for the locking stretch of the polymer chain networks, while α i define the additional coefficients (given in Table 1). The effective strain energy of homogenized composite depends on the invariants and interface defects volume fraction. It could further be divided into three independent contributions-of the matrix, of the interphase and of the particle, independently. There holds

Neo-Hookean
, I 1 -the first invariant of the strain tensor, µ-shear modulus Mooney-Rivlin W = C 1 I 1 − 3 + C 2 I 2 − 3 , C 1 , C 2 -empirical material constants, I 1 , I 2 -the first and the second invariants of the left Cauchy-Green deformation tensor It is aimed for composites with much softer matrix than the reinforcement. In such composites, deformation of stiff reinforcement is very small so that U p ξ ≪ U m ξ , U int ξ and the major contribution to deformation energy comes from hyper-elastic components. The relation between the stiffness tensor and the weakening coefficient w is determined with two steps (a) material coefficients are calculated on the basis of the FEM experiments solved sequentially with varying w and (b) these coefficients are replaced by a continuous polynomial with respect to w proposed as This polynomial is optimized with use of the Least Squares by maximization of correlation between the lab results and the computations as well as minimization of variance and error of the LSM. Such an approximation ensures continuous relation to w and ε 11 and introduces a very small error. The second formulation is based on the following approximation of the strain energy as a bivariate polynomial so that U e f f ,LSM = P(ε 11 , w) = n i=0 n j=0 a n,ij w·ε 11 , where polynomial coefficients are denoted by a ij and its rank by n. This approximation is next optimized in a two-step procedure. Firstly, a regression is applied for determination of polynomial coefficients a ξn with ranks between 1 and 12. Secondly, the optimum rank is selected that ensures a minimum LSM error. However, the principle objective of this study is to analyze numerically the deformation energy U for the hyper-elastic composite. Let X denote the position vector related to the non-deformed body configuration, x position vector related to the deformed configuration and u displacement vector. Let us define in the three-dimensional physical space for the given material the Green-Lagrange strain tensor as Let us define introduce the second Piola-Kirchhoff stress tensor as Energies 2020, 13, 2011 where J is the Jacobian of transformation from non-deformed to a deformed coordinate system. By definition of the hyper-elastic material we have the following relation: where W is an elastic potential dependent on the invariants of the strain tensor. One can derive the general constitutive relation for hyper-elastic isotropic material from the formula (15) as [30] where φ i is a derivative of the elastic potential with respect to the i-th invariant of the strain tensor Finally, we calculate elastic strain energy in a traditional way as It is done by assuming various forms of the elastic potential W corresponding to the constitutive models of hyper-elastic materials contrasted in numerical analysis and summarized in the table below. As it is seen, Arruda-Boyce theory adopted to the stochastic interface theory, includes higher powers of the strain tensor first invariant. So that, a precise initial FEM error analysis is extremely important to avoid any computational discrepancies during deformation energy determination; one could consider an application of some mesh adaptation procedure.
The Finite Element Method discretization of the entire RVE has been prepared with the use of 20-node three-dimensional iso-parametric finite element in three-dimensional physical space called C3D20 in ABAQUS. Let us introduce the global coordinate system x 1 , x 2 , x 3 and the local normalized coordinate system ξ 1 , ξ 2 , ξ 3 . The global coordinates of the i-th node after deformation we denote as where N k are the applied shape functions: for corner nodes: for side nodes: The coordinates of the element nodes in the local system ξ 1 k , ξ 2 k , ξ 3 k are shown in Figure 3.  We consider an iso-parametric finite element, so we can use the same shape function functions (20) to describe displacements increments where q is a vector containing the increments of the given degrees of freedom. Equation (24) can be further used to discretize the strain tensor increments ε defined previously in formula (13) and one arrives at The matrix B includes all spatial derivatives of the shape functions N. This allows for a determination of the stress tensor increments in the following way: The right hand side is expanded using random Taylor series of the following form: where   0, 11 ; eff Uw  is the mean value of effective energy under consideration, which means effective energy calculated for the mean value of the parameter w, and where odd order terms simply vanish. Such a procedure for the first four probabilistic moments of energy and for the given polynomial basis has been entirely programmed in MAPLE 2019, whose results are presented below. We consider an iso-parametric finite element, so we can use the same shape function functions (20) to describe displacements increments where ∆q is a vector containing the increments of the given degrees of freedom. Equation (24) can be further used to discretize the strain tensor increments ∆ε defined previously in formula (13) and one arrives at The matrix B includes all spatial derivatives of the shape functions N. This allows for a determination of the stress tensor increments in the following way: where C(ε) is strain-dependent constitutive tensor. A parametrization of this FEM model with respect to the interphase parameter w necessary in further Stochastic Finite Element Method implementation proceeds by a selection of the set of real values of this parameter taken from a certain neighborhood of its mean value, sequential recalculation of several FEM models with varying w and the Least Squares Method approximation of approximating polynomials independently for strains, stresses and energies (similarly to the formula (11)). Final calculation of the stochastic moments of the energy follows the classical integral definition of the kth central probabilistic moment, which is The right hand side is expanded using random Taylor series of the following form: where U 0,e f f (w; ε 11 ) is the mean value of effective energy under consideration, which means effective energy calculated for the mean value of the parameter w, and where odd order terms simply vanish. Such a procedure for the first four probabilistic moments of energy and for the given polynomial basis has been entirely programmed in MAPLE 2019, whose results are presented below.

Deterministic Numerical Experiments
Numerical experiments are done according to the homogenization method with use of a hexagonal Representative Volume Element (RVE). This RVE has unit dimensions and consists of three phases, i.e., the polymeric matrix that occupies 90% of the RVE volume, as well as a spherical particle and the surrounding interphase, that occupy remaining 10% of the RVE, 5% each. Linear elastic particles of carbon black C 60 have initial mechanical properties equal to E p = 10GPa and µ p = 0.3. The matrix and the interphase are modeled after the laboratory tests of Laripur LPR 5020. Fitting is provided with the Arruda-Boyce hyper-elastic potential and a set of parameters of C m 1 , λ m A , C int 1 , λ int A chosen from 10 different hyper-elastic models because it ensured optimal fitting to laboratory data. A surface-based tie constraint is selected for modeling a perfectly rigid contact. The interface defects volumetric ratio takes the value 0, when interphase is perfect and 1, when this interphase is composed from the defects only. Deformation energy of uniaxial stretch is computed with an implicit solver available in the FEM system ABAQUS. It is recovered from a series of cell problems, in which particle has been discretized using the hexahedral 20-noded C3D20 finite elements and the remaining phases-with the hybrid linear pressure C3D20H finite elements; the entire mesh shown in Figure 4 consists of about 150,000 finite elements. The full Newton algorithm has been applied in incremental analysis and the direct sparse multi-front equation solver has been used. The displacement of outer edges of the RVE is linearly increased in each increment of analysis. Size of increments ranges between 10 −5 and 5·10 −3 , while the outputs are recovered exactly 55 times during the deformation process for ε 11 ∈ (0.00, 0.275) for each FEM realization. The model of the matrix is optimized to best represent the mean value of the uniaxial tension laboratory tests of 10 specimens. The model used in further study is chosen by minimization of the total LSM error among 10 different hyper-elastic potentials. The fitting of material parameters is done by the Least Squares Method (LSM) with equal weights.

Deterministic Numerical Experiments
Numerical experiments are done according to the homogenization method with use of a hexagonal Representative Volume Element (RVE). This RVE has unit dimensions and consists of three phases, i.e., the polymeric matrix that occupies 90% of the RVE volume, as well as a spherical particle and the surrounding interphase, that occupy remaining 10% of the RVE, 5% each. Linear elastic particles of carbon black 60 C have initial mechanical properties equal to = 10 and = 0.3. The matrix and the interphase are modeled after the laboratory tests of Laripur LPR 5020. Fitting is provided with the Arruda-Boyce hyper-elastic potential and a set of parameters of ( 1 , , 1 , ) chosen from 10 different hyper-elastic models because it ensured optimal fitting to laboratory data. A surface-based tie constraint is selected for modeling a perfectly rigid contact. The interface defects volumetric ratio takes the value 0, when interphase is perfect and 1, when this interphase is composed from the defects only. Deformation energy of uniaxial stretch is computed with an implicit solver available in the FEM system ABAQUS. It is recovered from a series of cell problems, in which particle has been discretized using the hexahedral 20-noded C3D20 finite elements and the remaining phases-with the hybrid linear pressure C3D20H finite elements; the entire mesh shown in Figure 4 consists of about 150,000 finite elements. The full Newton algorithm has been applied in incremental analysis and the direct sparse multi-front equation solver has been used. The displacement of outer edges of the RVE is linearly increased in each increment of analysis. Size of increments ranges between 10 −5 and 5 • 10 −3 , while the outputs are recovered exactly 55 times during the deformation process for 11 ∈ (0.00, 0.275) for each FEM realization. The model of the matrix is optimized to best represent the mean value of the uniaxial tension laboratory tests of 10 specimens. The model used in further study is chosen by minimization of the total LSM error among 10 different hyper-elastic potentials. The fitting of material parameters is done by the Least Squares Method (LSM) with equal weights. Fittings and results of laboratory tests (marked by the green points) are presented in Figure 5. They include 10 different potentials, namely the reduced (2nd order) and full polynomials (2nd, 4th, 5th and 6th order) as well as Neo-Hookean [31,32], Arruda-Boyce [33], Yeoh [34] and Mooney-Rivlin [35,36] models. Arruda-Boyce potential has been selected for further analysis with the following parameters: 1, = 3.749 • 10 7 , = 7.00. It overestimates a little bit an overall value of the longitudinal stress at the strain level close to 11 = 0.275, but ensures at the same time an efficient overall estimation. As the confirmation one may notice that maximum stress coming from laboratory tests equals to Fittings and results of laboratory tests (marked by the green points) are presented in Figure 5. They include 10 different potentials, namely the reduced (2nd order) and full polynomials (2nd, 4th, 5th and 6th order) as well as Neo-Hookean [31,32], Arruda-Boyce [33], Yeoh [34] and Mooney-Rivlin [35,36] models. Arruda-Boyce potential has been selected for further analysis with the following parameters: C m 1,AB = 3.749·10 7 , λ m A = 7.00. It overestimates a little bit an overall value of the longitudinal stress at the strain level close to ε 11 = 0.275, but ensures at the same time an efficient overall estimation. As the confirmation one may notice that maximum stress coming from laboratory tests equals to E[σ m (ε 11 )] = 23.12 MPa and the one calculated from the Arruda-Boyce is equal   and for an increasing level of strain 11 . The direct output of the FEM simulations is marked as 'FEM', the functional approximation as 'LSM' and a bivariate polynomial approximation as 'BPO' on this figure. It illustrates that the deformation energy exponentially increases together with an additional increase of the strain level and it has a converse relation with for all levels of the strain. Effective deformation energy differs approximately by 40% between = 0.1 and = 0.9 for the highest considered strain level and this decrease gets smaller together with decrease of 11 .  Figure 6 reports a fluctuation of the effective deformation energy w.r.t. the interface defects volume fraction w and for an increasing level of strain ε 11 . The direct output of the FEM simulations is marked as 'FEM', the functional approximation as 'LSM' and a bivariate polynomial approximation as 'BPO' on this figure. It illustrates that the deformation energy exponentially increases together with an additional increase of the strain level and it has a converse relation with w for all levels of the strain. Effective deformation energy U e f f differs approximately by 40% between w = 0.1 and w = 0.9 for the highest considered strain level and this decrease gets smaller together with decrease of ε 11 .
Further, a relative error of the LSM and BPO is illustrated in Figure 7, which has been calculated as follows Err U e f f ,Approx = m i=1 U ξ,exact s,i (ε 11,i , w i ) − P n (ε 11 = ε 11,i , w = w i ) (29) Thanks to the series of previous computer experiments with this composite [37], an overall numerical error has been limited to the values equal or smaller than 1%, see Figure 7. It is a little bit expected that this error variations resulting from strain level fluctuations form quite regular patterns and typically decreases together with an increase of the strain level. An influence of the interface defects is of the accidental nature, where the extreme values are obtained for the defects of the smallest size and this deserves more detailed FEM tests, possibly with multiscale direct modelling of these defects. Such an analysis could be provided using the NURBS technique implemented together with stochastic perturbation-based formulation of the FEM [38,39]. Quite analogous observation concerns a coincidence of the BPO and the FEM series-they coincide almost perfectly while varying the strain level, whereas uniform changes of the parameter w results in accidental coincidence. Further, a relative error of the LSM and BPO is illustrated in Figure 7, which has been calculated as follows Thanks to the series of previous computer experiments with this composite [37], an overall numerical error has been limited to the values equal or smaller than 1%, see Figure 7. It is a little bit expected that this error variations resulting from strain level fluctuations form quite regular patterns and typically decreases together with an increase of the strain level. An influence of the interface defects is of the accidental nature, where the extreme values are obtained for the defects of the smallest size and this deserves more detailed FEM tests, possibly with multiscale direct modelling of these defects. Such an analysis could be provided using the NURBS technique implemented together with stochastic perturbation-based formulation of the FEM [38,39]. Quite analogous observation concerns a coincidence of the BPO and the FEM series-they coincide almost perfectly while varying the strain level, whereas uniform changes of the parameter w results in accidental coincidence.

Stochastic Numerical Analysis
The final study includes computation of stochastic characteristics of the composite effective deformation energy relative to statistical variability of the input random parameter . These analyses are performed thanks to the common usage of the systems ABAQUS and MAPLE 2019. They are performed for fluctuating expectation of the interface defects volume fraction in range of ( ) ∈ [0.1, 0.3, 0.5], for its increasing coefficient of variation ( ) ∈ (0, 0.25) and for different levels of the uniaxial strain within 11 ∈ (0, 0.275). Approximating optimal polynomials have 5 th or 6 th order, whereas some numerical studies available in the literature provide linear approximation for the same phenomenon [40]. The probabilistic characteristics analyzed in this study include the expected value

Stochastic Numerical Analysis
The final study includes computation of stochastic characteristics of the composite effective deformation energy U e f f relative to statistical variability of the input random parameter w. These analyses are performed thanks to the common usage of the systems ABAQUS and MAPLE 2019.
They are performed for fluctuating expectation of the interface defects volume fraction in range of E(w) ∈ [0.1, 0.3, 0.5], for its increasing coefficient of variation α(w) ∈ (0, 0.25) and for different levels of the uniaxial strain within ε 11 ∈ (0, 0.275). Approximating optimal polynomials have 5th or 6th order, whereas some numerical studies available in the literature provide linear approximation for the same phenomenon [40]. The probabilistic characteristics analyzed in this study include the expected value (Figure 8), coefficient of random dispersion (Figure 9), skewness ( Figure 10) and kurtosis (Figure 11) of the resulting U e f f (w, ε 11 ); skewness is output distribution asymmetry measure and kurtosis is output distribution concentration measure. These characteristics are determined separately with two independent methods, i.e., the fourth order iterative stochastic perturbation technique (SPT) and Monte-Carlo simulation (MCS) computed for 250,000 trials for any discrete point given in these figures. Computation is based on the full bivariate polynomial. Its coefficients are optimized with use of WLSM of equal weights and optimum rank is chosen among first ten orders. Equal weighing is preferred over Dirac weighting scheme because all the computations are based on a single bivariate polynomial and the expected value of w differs among realizations. It must be highlighted that computer resources required by the perturbation-based approach are much smaller than the ones necessary for Monte-Carlo simulation and time efficiency is also significantly better; the latter method is able to provide only discrete results that are marked on subsequent graphs with an asterisk.        Figure 11. Kurtosis of the effective deformation energy ( ). Figure 11. Kurtosis of the effective deformation energy κ U e f f . The expectations of the effective deformation energy E U e f f presented in Figure 8 are practically independent of the input random dispersion α(w), decrease together with an increase of expected value of defects volume fraction E(w) and increase together with an additional increase of ε 11 . A limited influence of α(w) on E U e f f is quite common for the Gaussian input parameters and was also recognized for the elastic case in [17]. An extent of decrease of E U e f f between extreme E(w) is comparable to the deterministic case presented on Figure 6.
Coefficients of variation α U e f f presented in Figure 9 are practically independent of ε 11 and increase together with an increase of the coefficient of variation of an input random variable α(w) as well as its expectation. Quite importantly, α U e f f for E(w) = 0.1 is nearly 10 times lower than for E(w) = 0.5, which means that the volume fraction of defects in the composite plays a crucial role not only for the energy itself but also its dispersion; this dispersion is smaller than the one of w but considerable knowing the volume fraction of defects (w = 0.1 means that V d corresponds to approx. 0.005 V). Correspondence of the two probabilistic methods for these basic probabilistic characteristics (E U e f f , α U e f f ) is perfect for the entire analyzed statistical scattering. This means remarkable time savings in computational analysis since a time effort of the stochastic perturbation technique is closer to a several solutions of the FEM problem than to the statistical large size repetitive solution of the homogenization problem typical for the Monte-Carlo analysis.
Higher order probabilistic characteristics, i.e., skewness β U e f f ( Figure 10) and kurtosis κ U e f f ( Figure 11) are both increasing their magnitude together with an increase of α(w) as well as E(w) and generally decreasing in magnitude together with an increase of ε 11 . Skewness remains always negative, whereas kurtosis-almost always positive. The level of the strain is important especially for low strains, for which the magnitudes of characteristics and their rate of change are the highest. Results of alternative probabilistic methods are still very good. Some differences are noticeable only for skewness with strains close to 0 and very high α(w) > 0.22 as well as for kurtosis with E(w) = 0.5 and α(w) > 0. 16. The resulting homogenized deformation energy is definitely not Gaussian because of a non-zero skewness and kurtosis-its PDF changes character together with a change of the level of the strain and this result is similar to other stochastic interphase models [41]. Characteristics are most sensitive to strain changes close to zero and are also not unique for E(w). An increase of E(w) always decreases expectation and increases dispersion as well as magnitudes of higher order characteristics. This means that interface defects are not only decreasing overall deformation energy of the composite but also increase a band of possible deformation curves. Moreover, higher order probabilistic characteristics computed using stochastic perturbation technique, even in the iterative mode [42], do not coincide with the statistical estimators of skewness and kurtosis computed via the Monte-Carlo scheme. They diverge from each other remarkably especially while increasing input coefficient of variation of the interphase parameter w.

Conclusions
It has been demonstrated here by numerical simulation that stochastic interface defects remarkably affect deformation process of the given composite. They considerably decrease a deformation energy of the homogenized composite U e f f at all the strain levels (up to 40%) and also cause its noticeable dispersion up to 5%. The resulting deformation energy cannot have Gaussian distribution and its PDF varies for different levels of the strain, especially close to 0. This means that application of a constant uncertain parameter independent of strain level for this composite is incorrect. The expected interface defects volume fraction E(w) has an influence on magnitudes of all stochastic characteristics of U e f f . The lower ones, i.e., expected value and coefficient of variation are increased and the higher ones, i.e., skewness and kurtosis are decreased together with its increase. The E U e f f is almost independent of α(w) and all higher probabilistic characteristics of U e f f always increase their magnitude when dispersion of the defects volume fraction α(w) increases. It is verified here that the proposed defective composite model is efficient in catching the degradation of composite structural response with an increasing volume fraction of interface defects; numerical error of this FEM approximation is always kept below a single percent. Further, it is proved here that the iterative stochastic-perturbation technique implemented in conjunction with the Finite Element Method is efficient and reliable for calculation of the first four probabilistic characteristics of deformation process in hyper-elastic particle-reinforced composites with uncertain interface defects. The coincidence of the proposed stochastic perturbation-based method with statistical analysis using classical Monte-Carlo simulation is almost perfect for the expected values and coefficients of variation and very good for the skewness and kurtosis. This is valid for expected value of interface defects volume fraction E(w) between 0.1 and 0.5, for its coefficient of variation α(w) as high as 0.25 and irrespective to the level of the strain ε 11 . Some divergence is visible only for kurtosis with E(w) = 0.5 and α(w) > 0. 16. This proves that the iterative stochastic-perturbation technique could replace a highly time and computational resource expensive Monte-Carlo technique well suited for discrete calculations and not especially for probabilistic processes. Further research investigations should focus on application of the stochastic techniques for truncated Gaussian or even non-Gaussian variables or processes to model stochastic effective stiffness of various nonlinear composites including interface defects.