Surrogate modeling of the effective elastic properties of spherical particle-reinforced composite materials

This paper focuses on the development of a surrogate model to predict the macroscopic elastic properties of polymer composites doped with spherical particles. To this aim, a polynomial chaos expansion based Kriging metamodeling technique has been developed. The training experimental design is constructed through a dataset of numerical representative volume elements (RVEs) considering randomly dispersed spherical particles. The RVEs are discretized using finite elements, and the effective elastic properties are obtained by implementing periodic boundary conditions. Parametric analyses are reported to assess the convergence of the scale of the RVE and the mesh density. The accuracy of the proposed metamodelling approach to bypass the computationally expensive numerical homogenization has been evaluated through different metrics. Overall, the presented results evidence the efficiency of the proposed surrogate modelling, enabling the implementation of computationally intensive techniques such as material optimization.


Introduction
The field of composite materials design has represented an area of rapid and steady growth since the early 1960s, with groundbreaking impacts across many scientific disciplines and industries [1][2][3]. Composite materials are defined as the artificial combination of two or more materials, usually known as phases, comprising from classical reinforced concrete or masonry to next-generation nano-modified composites. The simplest composite configuration combines two materials, a doping filler (inclusions) and the host matrix material. The theoretical estimation of the effective properties of composites is of critical importance to enable their design, minimizing the need for physical prototypes and experimentation and thus reducing costs and production times. However, the development of high-fidelity models of composite microstructures may be a challenging task involving extraordinary computational burden.
When designing high-strength composite materials, it is of pivotal importance to estimate the overall elastic constitutive tensor of the composite C * . To do so, a RVE is usually defined as the material domain containing a sufficient number of inclusions to statistically represent the composite as a whole. The dimensions of the RVE must be much smaller than the characteristic length over which the macroscopic loading varies in space. From a mathematical standpoint, let us consider two-phase RVE occupying a domain O ⊂ R 3 with boundary = ∂O. The elastic properties of the composite are defined by the geometric arrangement of the matrix O m and the inclusions O f , satisfying O m ∪O f = O with respective (fourth order) elastic tensors C m and C f . Let us denote ζ m and ζ f the volume fractions of the matrix phase and the filler, respectively, given by ζ s = |O s |/|O|, s = m, f , where |·| relates the volume of the constituent phases. The elastic response of the RVE to a certain boundary force field φ 0 and a volume force field g is governed by an elliptic steady-state problem [1,[4][5][6][7][8][9][10]: where φ is the displacement field, m ∪ f = with m ∩ f = ∅, and ν the outer unit normal to . The strain tensor ε is defined through the displacement field φ as: and it is related to the stress tensor by the generalized Hooke's law σ = Cε. The stiffness tensor of a material point x in the RVE is given by C as: Extracting the exact solution of the displacement field from Eq. (1) is a challenging task, and closed-form solutions can be only found for simplified or ideal composite microstructures. Nonetheless, in practical applications, it suffices to determine the overall behaviour of the composite at the macroscale, that is the asymptotic behaviour (homogenization) of Eq. (1) [11]. This is achieved through the homogenised version of Eq. (1) including the homogenised elastic tensor of the composite material C * as: Numerical homogenisation techniques such as the Finite Element Method (FEM) [12][13][14] constitute a popular approach to solve Eq. (3). In general terms, these techniques discretize the RVE into discrete elements, and solve the weak form of Eq. (3). Although these techniques enable a high-fidelity representation of the geometry of the microstructure, the computational time and computer resources may be substantial. This is particularly limiting when dealing with random arrangements of inclusions, which require the use of considerably dense discretization meshes. Such high computational demands may hinder or preclude their implementation in material optimization or uncertainty propagation analyses, which often require an elevated number of model evaluations.
In view of these limitations, we propose in this work the construction of a surrogate model to bypass computationally dense FEM models for the homogenization of composite materials. Specifically, we propose the use of polynomial-chaos expansion based Kriging metamodelling. The surrogate model is trained using a training dataset or experimental design constructed by a discrete set of evaluations of the numerical homogenization model covering the design space of the parameters of interest [15,16]. Once trained, thanks to the minimal computational cost involved in the evaluation of the metamodel, uncertainty propagation analyses are conducted to characterize the elastic properties of the composite in statistical terms. In this work, we focus on the analysis of epoxy composites doped with glass microsphere particles. These composites exhibit remarkable mechanical properties such as high bending/compression strength, bulk modulus, wear resistance and coefficient of friction, finding a wide range of applications in numerous fields such as automotive, aeronautics and biotechnology [17][18][19].
The remainder of this paper is organized as follows. Section 2 outlines the FEM homogenisation approach used to estimate the effective properties of RVEs of epoxy doped with glass microsphere particles. Section 3 shows the theoretical formulation of the proposed surrogate model. Section 4 presents the numerical results and discussion and, finally, Sect. 5 concludes this work.

Numerical homogenisation of glass microsphere/epoxy resin composites
Several strategies can be followed for the numerical homogenisation of composite materials. The simplest one consists in the use of periodic unit cells [20], containing a canonical arrangement of particles which are assumed to replicate periodically throughout the composite material. The mesh density required to discretize such RVEs is limited and so is the computational burden involved in the homogenisation. Nevertheless, most composite materials present a certain degree of randomness in the dispersion of the reinforcing fillers. Although there are 2D models capable of representing this randomness [21], it is often necessary to implement more computationally intense 3D random RVEs [22]. In this case, different simplifying assumptions can be also adopted, including RVEs allowing particles to intersect with the RVE boundaries (see e.g. [23]), and the more general case of intersecting particles [24]. In this work, we implement a general cubic 3D model with allows particle intersection with the cell boundaries. For the construction and analysis of the material microstructure, a combination of scripts generated in MATLAB environment and in the commercial FEM code ANSYS are used. Specifically, the adopted methodology for defining the geometry of the cubic RVE involves the following four steps: 1. Given the particle radius, generate a particle center coordinates (x, y, and z). That center is allowed to fall outside of the volume. 2. Check if the particle cuts any of cubic surfaces. If so, symmetric particles-one, three or seven, depending if the original one cuts one, two or three cell faces, respectively-are generated in order to impose RVE periodicity. 3. Check if all particles generated in the previous step intersect the previous existing ones. If this is not the case, the proposal particle is accepted, otherwise return to step 1. 4. Repeat until the desired filler volume fraction is achieved.
Once the geometry is set up, it is discretized in ANSYS using 4-nodes linear thetraedral solid elements (SOLID 285). A sample of one of the generated RVEs is shown in Fig. 1a and the corresponding mesh is depicted in Fig. 1b.
To compute the effective properties of composite, it is assumed that the material can be defined as the periodic replication of the previously defined RVE as shown in Fig. 2a. Therefore, to determine the effective elastic properties of one single isolated RVE, periodic boundary conditions must be applied. The general periodic boundary conditions on the cell faces of a RVE are given by [25]: whereε i j denote the volume average strains, and v i represents the local periodic part of the displacement components u i on the boundary surfaces. The latter displacement components are generally unknown and depend upon the applied loading. Indices i and j denote the global Cartesian directions. In the case of square RVEs like the ones used in this work and sketched in Fig. 2a; Eq. (4) takes a more explicit expression. Fig. 2a. Then, the displacements on a pair of opposite boundary cells (with their normal along the x j axis) read:

Consider the notation of the cell surfaces
where indexes K + and K − indicate the displacements along the positive and negative x j direction, respectively. Local fluctuations v K − i and v K + i must be identical on every two opposing faces to comply with the periodicity of the RVE. Therefore, the local displacement components can be dropped from the formulation by the difference between the expressions in Eq. (5), leading to: Therefore, the RVE can be subjected to a desired strain state by imposing proper displacements on its boundary surfaces. Then, the volume average stressesσ i j and strainsε i j in the RVE can be computed as: with V being the volume of the RVE. Then, the i j-th component of the elastic tensor can be directly estimated as C * i j =σ i j /ε i j . Since, in general, the elastic tensor of an anisotropic material can be written following Voigt's notation as a symmetric 6 × 6 matrix as [26]: a total of six virtual tests are required to fully characterize the homogenised tensor C * (3 pure dilations and 3 pure distortions). Alternatively, we can rewrite Eq. (9) in terms of the compliance matrix S * as: ⎛ so the engineering constants of interest are given in terms of S i j [26]: In the particular case of isotropic materials (like is is the case of a RVE made of an isotropic matrix material doped with a sufficient number of isotropic spherical particles), the components C * i j in Eq. (9) must satisfy: where C * 12 = λ * and C * 44 = μ * are the usual Lamé moduli. Note that he Lamé constants are related to the Young's modulus E * and Poisson's ratio ν * of the homogenised material as follows [27]: The boundary conditions in Eq. (6) have been implemented in ANSYS by coupling the displacements of opposite nodes and opposite boundaries. For each pair of displacement components at two corresponding nodes with identical in-plane coordinates (within a certain tolerance) on two opposite cell faces, the corresponding constraint condition according to Eq. (6) is imposed. As an example, Fig. 2b shows the constraint equations to be defined for a pair of nodes located on the opposite surfaces A − and A + . Interested readers may find further information on the computational implementation of the boundary conditions in Ref. [25].

Polynomial-chaos expansion based Kriging
Consider a mathematical model M which maps a M-dimensional input parameter space to a 1-dimensional output space.
Due to uncertainties in the input vector, x is considered a random variable with known probability distribution. A surrogate model or metamodel is a mathematical function M which aims to emulate the original model M , but at a lower computational cost. The general procedure to construct a surrogate model can be summarized in the following steps: 1. Sample a training set T and a validation set V covering the parameter design space: In this work, we propose the combination of Kriging and polynomial chaos expansion. Kriging metamodelling has been widely used in the literature to bypasss models exhibiting important local effects and/or nonlinearities, while polynomialchaos expansion has proved highly efficient to represent the global response of models with no marked local effects. In this work, we employ a combination of Kiring and polynomial chaos expansion as a computationally efficient metamodel with both global and local modelling capabilities.
• ψ (i) α i are orthogonal polynomials depending on the stochastic nature of the input variables x (e.g. Legendre or Hermite polynomials). In this work Legendre polynomials will be used, as variables under study are supossed to follow an uniform distribution.
Although the expression (15) can be proved exact for an infinite number of polynomials, in practice only finite terms of α∈N M a α α (x) can be computed. As a consequence, different strategies can be taken into account to truncate the polynomial series. The simplest one consists in selecting all the polynomials whose total degree |α| = M i=1 α i belongs to the set where the cardinality of A M, p is equal to Nonetheless, when M and p are big enough, this procedure of polynomial selection may lead to computing a large number of coefficients, resulting in large computational burdens. As an alternative, an hyperbolic truncation scheme can be employed. This approach consists in selecting all multi-indices with q-norm α q = M i=1 α q i 1 q less or equal to p, i.e.: This strategy has been proved efficient due to the fact that high interaction terms tend to have coefficients close to zero in many practical applications. Note that, when q = 1, this scheme corresponds with the previous Although substantial cost reductions can be achieved using the hyperbolic truncation scheme, the number of coefficients in the expansion may still be considerable. A more cost-efficient solution can be obtained by using the adaptive least-angle regression (LAR) algorithm [29]. LAR constructs a set of expansions incorporating an increasing number of basis polynomials α , from 1 to P = card A M,p,q . Then, the resulting sequence of index sets is used to construct different expansions and the best metamodel is selected by a cross validation procedure. Finally, the expansion coefficients a = {a α , α ∈ A M, p ⊂ N M } are obtained by minimizing the expectation of the least squares errors: where x (i) ∈ T , i = 1, . . . , N . Once Eq. (16) has been solved, the resulting PCE surrogate model can be written as follows:

Kriging
The Kriging model assumes that the response of a computational model M (x) can be modelled as the sum of a deterministic trend and a stochastic process as: where T (x) is a regression model, often termed the trend of the Krigng model, and Z (x) is a zero-mean stochastic process. An interpretation of Eq. (18) is that deviations from the regression model may resemble a sample path of a properly chosen stochastic process [30]. Depending on the nature of the trend part T (x), we can distinguish simple, ordinary on universal Kriging models. Simple Kriging assumes that T (x) is a known constant value, Ordinary Kriging supposes that T (x) is an unknown constant value and, finally, Universal Kriging considers that T (x) can be described by a general linear regression model: Ordinary Kriging is one of the most popular metamodelling approaches in the literature [31], although a proper definition of the trend part in Universal Kriging can lead to more accurate results, as evidenced in the numerical results presented hereafter. Note that simple and ordinary Kriging are special cases of Universal Kriging.
Like in the case of PCE, calibrating a Kriging surrogate model requires solving an optimization problem. The stochastic process Z (x) is determined by an autocorrelation function: depending on the distance between points x and x and on some hyperparameters θ ∈ R M to be determined. This function describes the correlation between points of the input space R M . In general, the correlation decreases with distance |x − x |, and larger values of θ leads to a faster decreasing. In this work, Gaussian auto-correlation functions are considered as [32]: The computation of hyper-parameters θ is usually performed by either the maximum-likelihood-estimation [33] (ML) or by the leave-one-out cross validation [34] (CV), which are described by the following minimization problems: It is no clear which approach leads to better results [34], depending the answer on the specific problem under study. In this work, both techniques have been implemented achieving very similar results. In order to solve the optimization problems in Eqs. (22) or (23), different techniques can be employed. In general, optimization techniques can be classified as local (e.g. gradient-based) or global optimization techniques (e.g. genetic or differential evolution algorithms). While local optimisation algorithms are generally less computationally demanding, there may be difficulties to find global extrema and the solutions often get stuck in local minimum. For this reason, a genetic algorithm has been used in this work to solve the problems above.

Polynomial chaos expansion based Kriging
Polynomial chaos expansion based Kriging is a surrogate model that combines in a simple way the two previous approaches, taking advantage of both techniques' characteristics. Following [35], a PCE-Kriging (PCK) surrogate model is an Universal Kriging model where the trend part T (x) is particularized as the truncated PC expansion using the LAR procedure as follows: Then, building the PCK metamodel consists in two steps. First, the optimal set of orthonormal polynomials α is determined using the LAR procedure (for α ∈ A the truncation set). Subsequently, the hyperparametersθ of the stochastic process are determined considering the previous PCE expansion as the trend part. As anticipated above, PCE approximates well the global behavior of the computational model, whereas Kriging manages the local variability of the model output [35].

Numerical results and discussion
In this section, numerical results obtained after the homogenization process are presented. Firstly, a preliminary convergence study to evaluate the required mesh density and the scale of the RVE is shown. Afterwards, the training of the PCE-Kriging metamodel is presented and numerical results are reported to evaluate its computational efficiency and accuracy. The material properties of epoxy and glass spheres have been taken from Ref. [36] and listed in Table 1.

Convergence analysis of the RVE length and mesh density
In order to assess the quality of the numerical RVE, two important aspects need to be carefully appraised: the mesh density and the size of the RVE. The discretization of the RVE must be fine enough to ensure that the isotropic homogenised constitutive tensor C * is mesh independent. On the other hand, the length of the RVE must be defined in such a way that it is effectively representative of the composite material as a whole. Therefore, the size of the RVE must be defined after a convergence analysis, ensuring that the estimates of the homogenised elastic tensor reach convergence values. Figure 3 reports the convergence analysis of the mesh density, carried out on a 100 μm edge cell containing 31 particles, showing the Young's moduli E 1 , E 2 and E 3 versus the total number of nodes in the mesh. Following the constitutive relation (10), the Young's modulus in the i-th (i = 1, 2, and 3.) direction is represented by the ratio between the stress and the strain along this direction, that is E i = σ i /ε i . Seven different meshes with increasing densities have been considering, including 6525, 12138, 42519, 90443, 158881, 231181 and 318727 nodes. Assuming the finest mesh attains the real values, the 4th one-150000 nodes-presents differences with the last one-320000 nodes-less than 0.5%.
A similar procedure has been carried out with the purpose of finding the minimal cell size that assures the isotropy of the RVE. To do so, four different cubic cells have been created, with edge lengths of 120, 180, 240 and 300 µm respectively. In this light, ratios E 1 /E 2 , E 1 /E 3 and E 2 /E 3 are furnished versus the total number of nodes in the RVE in Fig. 4. The composite material can be assumed as isotropic when these ratios are close to 1 (within a certain tolerance). Due to the 3D nature of the numerical model, increasing the edge sizes of the cell under study dramatically raises the computational burden in the homogenization. For instance, while the RVE with edge length 180 µm needs about 35 minutes to perform the homogenization, the RVEs with edge lengths of 240 µm and 300 µm require more than two and ten hours, respectively. As a consequence, a RVE size of 240 µm has been chosen as a trade-off between accuracy and computational cost. Note that even the smallest cells also show a good behaviour in terms of isotropy, although a greater size has been chosen in order to minimize numerical approximation errors.

Surrogate modeling
Following the theoretical framework introduced in Sect. 3.3, a PCE-Kriging metamodel has been constructed using 20 training points and validated with 10 independent  Fig. 4 Ratios of elastic moduli E i /E j , i, j = 1, 2, 3 versus the edge length of cubic RVEs of random glass microsphere/epoxy composite points. The variables included in the model are the particle size, ranging from 20 to 3 µm, and the volume fraction of particles, ranging from 0 to 20%. Both the training set T and the validation set V have been selected by means of the Latin Hypercube Sampling [37] procedure. To evaluate the accuracy of the surrogate model, the coefficient of determination R 2 and the normalized average absolute Error (NAAE) are considered as: where V = x (1) , . . . x (K ) is the validation set, and are the outputs of the numerical model and the surrogate model evaluated on V , respectively. The estimates of the surrogate model of the effective Young's modulus and Poisson's ratio random glass microsphere/epoxy composite are shown as functions of the particles' size and volume fraction in Fig. 5. The validation datapoints obtained by the numerical homogenisation of the RVE are also depicted with red scatter points. Overall, the Young's modulus of the homogenised material is determined with error measures of R 2 = 0.9991 and NAAE = 0.0237, thus evidencing the great performance of the PCK metamodel. It is noticeable in Fig. 5 that almost all the variability of the material properties is due to the volume fraction of the reinforcement particles, whit the particles' size has almost negligible effect. This agrees with some previously reported results in the literature (see e.g. [38,39]). With regard to the estimation of the effective Poisson's ratio, error measures of R 2 = 0.9983 and NAAE = 0.0320 are obtained, which also proves high accuracy in the metamodeling. It is also important to remark the high computational efficiency attained by the developed surrogate model. Specifically, while the numerical FEM takes about three hours for one single homogenisation, the metamodel only takes a few seconds to perform 30,000 model evaluations.

Concluding remarks and future work
In this study, highly accurate PCK surrogate models have been developed to bypass the computationally intensive numerical homogenization of the elastic properties of random glass microsphere/epoxy composites. To do so, numerical RVEs considering random arrangements of glass microspheres embedded in epoxy have been developed. The effective elastic properties of the RVEs have been obtained through numerical homogenisation and periodic boundary conditions. On this basis, an experimental design of 20 training points and 10 validation points have been obtained. Then, the training dataset has been used to train the PCK surrogate models considering the homogenised Young's modulus and Poisson's ratio as independent variables, and the particle's size and volume fraction as design parameters. The reported results evidence that, although a limited training dataset comprising only 20 training points has been used, the developed surrogate model exhibits high accuracy and extraordinarily low computational evaluation times. The reported results demonstrate that the volume fraction of the reinforcement particles critically determine the overall elastic properties of the composite.
Ongoing research is focusing on the incorporation of more design variables in the metamodeling, including the elastic properties of the constituent phases. Additionally, future work will address the consideration of randomness and uncertainty in the design variables. This will allow to evaluate the propagation of uncertainties in the definition of the constituent properties to the effective elastic behavior of random glass microsphere/epoxy composites.
In addition, the proposed scheme can be applied to the computationally intensive molecular dynamics (MD) simulation [40]. A typical atomistic MD system computes the forces between 100, 000 atoms, and through solving the Newton law of motion, updates their position and velocity. This process must be replicated millions of times and requires several days to complete on a standard HPC node. Communication algorithms adapted to MD have been evaluated in terms of efficiency by using d-meshes [41] or by an optimal dimension reduction scheme of the full original parameters space [42]. In view of this, the calculation potential of PCK applied to MD simulation constitutes a future work proposal, not only to characterize the macromolecular response but also to analyze the relative importance of the input random variables involved through an metamodel global sensitivity analysis study.
Another relevant chemical problem is the study of potential energy surfaces (PES). A PES is a mathematical relationship between the energy of a molecule (or a set of molecules) and its geometry. However, computing the energy of such system with sufficient accuracy implies excessive computational requirements. In order to bypass this problem, Kriging technique has recently been used to reduce the computational burden of such problems [43]. However, we believe that PCK approach may could be of interest in this type of problems to obtain higher quality surrogate models, taking also advantage of the good properties of both PCE and Kriging itself.