Stochastic Finite Element Analysis Framework for Modelling Mechanical Properties of Particulate Modified Polymer Composites

Polymers have become indispensable in many engineering applications because of their attractive properties, including low volumetric mass density and excellent resistance to corrosion. However, polymers typically lack in mechanical, thermal, and electrical properties that may be required for certain engineering applications. Therefore, researchers have been seeking to improve properties by modifying polymers with particulate fillers. In the research presented herein, a numerical modeling framework was employed that is capable of predicting the properties of binary or higher order composites with randomly distributed fillers in a polymer matrix. Specifically, mechanical properties, i.e., elastic modulus, Poisson’s ratio, and thermal expansion coefficient, were herein explored for the case of size-distributed spherical filler particles. The modeling framework, employing stochastic finite element analysis, reduces efforts for assessing material properties compared to experimental work, while increasing the level of accuracy compared to other available approaches, such as analytical methods. Results from the modeling framework are presented and contrasted with findings from experimental works available in the technical literature. Numerical predictions agree well with the non-linear trends observed in the experiments, i.e., elastic modulus predictions are within the experimental data scatter, while numerical data deviate from experimental Poisson’s ratio data for filler volume fractions greater than 0.15. The latter may be the result of morphology changes in specimens at higher filler volume fractions that do not comply with modelling assumptions.


Introduction
Having low volumetric mass density and high environmental degradation resistance (e.g., corrosion), while being low-cost, polymers are attractive for a variety of engineering applications [1][2][3]. Yet, limited mechanical, thermal, and electrical properties often limit the usability of polymers. This problem has motivated material designers to enhance polymer properties by adding appropriate types of fillers [4]. A variety of filler types and morphologies, in addition to the chosen filler volume fractions (FVF), are available for modifying polymers. The multitude of micro and nano-scale particles, and combinations thereof, forming not only single but binary, ternary, and even higher order modified polymers, create a need for modeling and experimental processes to guide the design of filler-modified polymers with well dispersed and distributed particles [5,6]. As such, the development of filler modified polymer composites is a costly and time-consuming process that is challenging to perform [7] by material designers.
Various studies have been reported for predicting the properties of advanced materials [8][9][10][11]. Eshelby pioneered the characterization of particulate composites by performing stress field analysis with an ellipsoidal inclusion [12]. The Halpin-Tsai empirical model, as shown in Equations (1) and (2), was proposed for predicting the modulus of elasticity for a variety of particle geometries [13].
where E f and E m are the modulus of elasticity of filler and matrix respectively, φ is the volume fraction of filler, and ζ and E c are correspondingly a shape factor and the effective modulus of elasticity. The approximation by Mori and Tanaka [14] is an extension of the Eshelby solution. Other analytical methods include modelling approaches by Benveniste [15], Christensen and Lo [16], and Torquato [17].
Some analytical approaches provide estimates for the upper and lower bounds of the modulus of elasticity with respect to the filler volume fraction. The Voight [18] and Reuss [19] models, as indicated by Equations (3) and (4), are the most basic upper and lower bounding techniques, respectively. The variational approach developed by Hashin and Shtrikman [20] is another popular analytical technique for predicting the bounds for the effective modulus of elasticity.
where E U and E L are the upper and lower bound effective modulus of elasticity, respectively. In addition to randomly oriented particles, modified polymers with aligned particles have gained interest among material designers for improving mechanical and thermal properties [21]. High strength filler modified polymers [22][23][24] and sensors [25,26] are example applications that require filler alignment to create explicit anisotropic mechanical properties. Carbon nanotubes, for example, are a popular filler type being considered for creating anisotropic materials for a variety of applications [27,28].
Analytical approaches are widely recognized as expedient for predicting mechanical properties of filler modified polymers. However, these methods often lack in accuracy, especially for high FVF. Experimental studies, on the other hand, provide direct information on material properties. However, being time-consuming and costly, experiments by themselves are an unsatisfactory solution for researchers designing advance materials [29]. In light of the aforementioned complexities associated with filler modified polymers, and a rapidly growing number of opportunities for these composites in the industry, it is imperative to identify alternative approaches that allow for the efficient and effective prediction of material properties [30][31][32].
Stochastic methods are prominently used in reliability analysis, which requires analyzing multiple input variables and predicting outcomes with a suitable level of accuracy. Stochastic methods have been powerful in other areas as well, such as complicated financial and forecasting models that involve several random and difficult to assess input variables. There are different methods of stochastic analyses, some relying on analytical approaches for predicting outcomes, e.g., first and second-order reliability methods. Using statistical principles, these methods enable predicting outputs from randomly determined input variables, which often follow certain types of random patterns or probability distributions. Stochastic methods thus consider input data variation or fluctuation and enable precise outcome estimation, even for complex systems. It is important to recognize that in some applications, uncertainty in input data may be dependent on each other. Such problems would be very difficult to predict, however stochastic analysis enables simulating the outcome by performing a large number of simulations known as stochastic projections [33].
It is convenient to consider explicit input parameters when solving engineering problems. As such, explicit outputs are also simulated. However, strictly speaking, explicit input values do not exist for material properties, boundary conditions, and structural geometries, since such data may be subject to variation or fluctuation. Therefore, considering input parameters that are explicit can lead to generating hidden errors in final solutions and outcomes that are far from reality. Alternatively, by considering uncertainty for input parameters, robust results can be obtained. Moreover, the number of input parameters that influence output values increases as engineering models become more comprehensive. Consequently, material designers are faced with multiple challenges, which highlights the importance of developing innovative simulation methods.
Stochastic finite element analysis (SFEA) is a statistical method developed recently for solving sophisticated engineering problems that are impossible to solve using analytical methods [34,35]. As the name implies, SFEA is capable of considering uncertainty for model input parameters. In the present context of filler-modified polymers, several random variables need to be considered for material design, e.g., particle locations, orientations, dimensions, and properties. All of these parameters are known to affect outcomes appreciably. Given the complex nature of this problem, SFEA is considered a viable approach for material properties prediction.
Two different options are pursued connecting stochastic analysis with finite element analysis (FEA) in order to perform SFEA on engineering models. The first option entails the use of dedicated commercial FEA software with stochastic analysis and scripting capability [36][37][38]. Stochastic analysis and scripting capabilities are hereby integrated into the software, making it convenient for a user to perform SFEA since the software applies uncertainty to input variables. However, this first option is limited in terms of flexibility to incorporate specific modeling requirements. The second and arguably more powerful option involves advanced scripting for connecting stochastic analysis with FEA, providing maximum modeling flexibility, e.g., in terms of applying uncertainty to input variables, choosing a suitable pseudorandom number generator and performing customized FEA processes, examples of which are demonstrated by reference [39]. However, the second option is typically more challenging for users to perform. In addition, with greater opportunities for capturing complexities in terms of model features and input data variation, the second option is prone to become computationally expensive, necessitating the use of high-performance computational resources. Recently, the present authors developed an SFEA framework for predicting the effective thermal conductivity of particulate polymer composites [40]. The interested reader is referred to this article for detailed information about the developed SFEA framework, which is briefly summarized in the subsequent section. In general, this approach is capable and intended of modeling filler particle sizes ranging from nano-to micro-scales. However, especially for nano-scale fillers, great care must be taken to properly identify and implement the particle-to-particle and particle-to-matrix interactions in order to correctly capture the effective material properties. For the present study, the SFEA framework was employed to explore mechanical properties, i.e., elastic modulus, Poisson's ratio, and coefficient of thermal expansion (CTE), for a composite with size-distributed spherical particles. Specifically, a composite consisting of an epoxy matrix with micro-sized spherical glass particles was studied and numerical results were contrasted with analytical predictions and experimental data from the technical literature.

Stochastic FEA Framework
The developed modeling framework employs Monte Carlo simulation techniques and computes outcomes using statistical analysis, which necessitates large numbers of model iterations. A customized stochastic analysis process in conjunction with parametric FEA was thus developed that automated the process of applying uncertainties to input variables by connecting several program modules using multiple programming languages. The illustration in Figure 1 is the high-level architecture of the algorithm used for the modeling framework. In the following, the various modules depicted in Figure 1, including associated acronyms, will be highlighted. routing information between the various modules used in the framework. A customized "Front End" form was created that enabled capturing of all input parameters required for the analysis and storing them in the database ('DBMS'). Input parameters captured by the "Front End" form included: Size information for the modeled representative volume element (RVE), FVF that were to be analyzed, filler and matrix properties, filler particle size distribution, boundary conditions, settings controlling the contact behavior for filler particles and matrix, and parameters for finite element (FE) mesh generation. An Open Database Connectivity (ODBC) method was used for developing the Database Management System (DBMS). The ODBC approach enabled creating an operating routine that was independent of the database and provides a high degree of flexibility for accessing the database at any time in the analysis. One task of the DBMS was the transfer of input parameters to the Monte Carlo Simulation (MCS) module in order to perform the SFEA.
Using VBA programming language, the MCS module was developed having a tabulated structure. The module enabled storing input parameters and results created by the model generation and analysis sub-processes. The algorithm as shown in Figure 2 is repeated for each specified FVF until statistical objectives are satisfied. The MCS module computes mean values, standard deviations, and variances, which were used as objective values for terminating the framework's iteration. Once acceptance criteria were satisfied, e.g., an explicit standard deviation or variance, the MCS stops iterating and calculates the final effective mechanical properties. Results calculated in each iteration and mean values are saved to the database, which can be accessed through the "Front End". Visual Basic for Applications (VBA) programming language was employed for connecting and routing information between the various modules used in the framework. A customized "Front End" form was created that enabled capturing of all input parameters required for the analysis and storing them in the database ('DBMS'). Input parameters captured by the "Front End" form included: Size information for the modeled representative volume element (RVE), FVF that were to be analyzed, filler and matrix properties, filler particle size distribution, boundary conditions, settings controlling the contact behavior for filler particles and matrix, and parameters for finite element (FE) mesh generation.
An Open Database Connectivity (ODBC) method was used for developing the Database Management System (DBMS). The ODBC approach enabled creating an operating routine that was independent of the database and provides a high degree of flexibility for accessing the database at any time in the analysis. One task of the DBMS was the transfer of input parameters to the Monte Carlo Simulation (MCS) module in order to perform the SFEA.
Using VBA programming language, the MCS module was developed having a tabulated structure. The module enabled storing input parameters and results created by the model generation and analysis sub-processes. The algorithm as shown in Figure 2 is repeated for each specified FVF until statistical objectives are satisfied. The MCS module computes mean values, standard deviations, and variances, which were used as objective values for terminating the framework's iteration. Once acceptance criteria were satisfied, e.g., an explicit standard deviation or variance, the MCS stops iterating and calculates the final effective mechanical properties. Results calculated in each iteration and mean values are saved to the database, which can be accessed through the "Front End". The Random Number Generator module (RNG, see Figure 1) was developed using the general mathematical programming environment MATLAB, enabling the generation of random numbers required for input variables. In the present study, these variables are Cartesian coordinates for particle locations (X, Y, Z) and the particle diameter. Particles were thus randomly distributed and dispersed within the RVE. The RNG module, illustrated in Figure 3, also performed collision detection between particles inserted into the RVE. A routine was defined that determined the distance for each particle newly inserted into the RVE with respect to all other particles already contained in the RVE. If no particle collision was detected, the new particle was accepted into the The Random Number Generator module (RNG, see Figure 1) was developed using the general mathematical programming environment MATLAB, enabling the generation of random numbers required for input variables. In the present study, these variables are Cartesian coordinates for particle locations (X, Y, Z) and the particle diameter. Particles were thus randomly distributed and dispersed within the RVE. The RNG module, illustrated in Figure 3, also performed collision detection between particles inserted into the RVE. A routine was defined that determined the distance for each particle newly inserted into the RVE with respect to all other particles already contained in the RVE. If no particle collision was detected, the new particle was accepted into the RVE; otherwise, the particle was rejected. This process continued until the specified FVF was satisfied, at which point the RNG module stops providing randomized particle data. detection between particles inserted into the RVE. A routine was defined that determined the distance for each particle newly inserted into the RVE with respect to all other particles already contained in the RVE. If no particle collision was detected, the new particle was accepted into the RVE; otherwise, the particle was rejected. This process continued until the specified FVF was satisfied, at which point the RNG module stops providing randomized particle data.
The composition of particulate filler materials may follow a certain particle size distribution. Simply using an average particle size may not reflect the behavior of the composite material. Hence, the SFEA framework was developed to account for particle size data that conforms to a prescribed size distribution using a data binning approach. The interested reader is referred to reference [40] for details on the algorithm that was included to conform particle dimensions to explicit size distribution, and for the influence on computational performance when adding size-ordered particles to the RVE. The latter was implemented in the SFEA framework. For the model generation and analysis sub-processes, a customized parametric FEA platform was developed using the commercial FEA software ANSYS Workbench (Version 19, ANSYS Inc., Canonsburg, PA, USA). This FEA platform as depicted in Figure 4 was realized in conjunction with The composition of particulate filler materials may follow a certain particle size distribution. Simply using an average particle size may not reflect the behavior of the composite material. Hence, the SFEA framework was developed to account for particle size data that conforms to a prescribed size distribution using a data binning approach. The interested reader is referred to reference [40] for details on the algorithm that was included to conform particle dimensions to explicit size distribution, and for the influence on computational performance when adding size-ordered particles to the RVE. The latter was implemented in the SFEA framework.
For the model generation and analysis sub-processes, a customized parametric FEA platform was developed using the commercial FEA software ANSYS Workbench (Version 19, ANSYS Inc., Canonsburg, PA, USA). This FEA platform as depicted in Figure 4 was realized in conjunction with scripting in IronPython programming language. Two distinct submodules were developed in order to create the parametric three-dimensional model geometry and subsequently the full FEA model in a two-step process. In the first step, ANSYS DesignModeler in conjunction with scripting in JAVA programming language was used to read the data generated by the RNG module and create geometries for particles and the RVE. The resulting geometric representation was then transferred to the FEA modeling environment. For the second step, the FEA environment was developed using ANSYS Mechanical, again in combination with JAVA scripting, which facilitated an automated process of reading data from the database and creating the FEA model, including the application of material properties, mesh generation, and extraction of FEA results. Subsequent to mesh generation, the FEA environment submodule performed a convergence study by refining the mesh and extracting results to check whether results convergence was satisfactory. Finally, results obtained with appropriately refined meshing were transferred to the MCS module and saved for further statistical analysis. process of reading data from the database and creating the FEA model, including the application of material properties, mesh generation, and extraction of FEA results. Subsequent to mesh generation, the FEA environment submodule performed a convergence study by refining the mesh and extracting results to check whether results convergence was satisfactory. Finally, results obtained with appropriately refined meshing were transferred to the MCS module and saved for further statistical analysis.

Model Details
The developed SFEA framework was herein employed to predict the mechanical properties of particulate polymer composites with randomly distributed spherical particles under static-structural condition. Specifically, spherical glass bead particles embedded in epoxy were considered. Properties were adopted from the technical literature accordingly [41,42], see Table 1. Materials were treated as linear elastic. The particle size distribution illustrated in Figure 5 was adapted from reference [43] assuming that all particles fall within a band ranging from 0 to 50 μm and a data binning approach with 5 bins was used to have particle sizes emulate the given distribution. The size of the cubical RVE was set at 400 μm based on preliminary size effect studies.

Model Details
The developed SFEA framework was herein employed to predict the mechanical properties of particulate polymer composites with randomly distributed spherical particles under static-structural condition. Specifically, spherical glass bead particles embedded in epoxy were considered. Properties were adopted from the technical literature accordingly [41,42], see Table 1. Materials were treated as linear elastic. The particle size distribution illustrated in Figure 5 was adapted from reference [43] assuming that all particles fall within a band ranging from 0 to 50 µm and a data binning approach with 5 bins was used to have particle sizes emulate the given distribution. The size of the cubical RVE was set at 400 µm based on preliminary size effect studies. Table 1. Properties of filler particles and polymer matrix.
Three-dimensional 10-node quadratic tetrahedral structural solid elements (SOLID187) were used for meshing both particles and matrix. Figure 6 illustrates an example mesh generated for the RVE and the particles included within it. The interface between the particles and matrix was defined by employing three-dimensional 8-node surface-to-surface contact elements (CONTA174 and TARGE170). In the present study, particle-matrix interfacial contact was assumed to constitute perfect bonding, and hence, no relative displacement at the interface was allowed. Three-dimensional 10-node quadratic tetrahedral structural solid elements (SOLID187) were used for meshing both particles and matrix. Figure 6 illustrates an example mesh generated for the RVE and the particles included within it. The interface between the particles and matrix was defined by employing three-dimensional 8-node surface-to-surface contact elements (CONTA174 and TARGE170). In the present study, particle-matrix interfacial contact was assumed to constitute perfect bonding, and hence, no relative displacement at the interface was allowed.
used for meshing both particles and matrix. Figure 6 illustrates an example mesh generated for the RVE and the particles included within it. The interface between the particles and matrix was defined by employing three-dimensional 8-node surface-to-surface contact elements (CONTA174 and TARGE170). In the present study, particle-matrix interfacial contact was assumed to constitute perfect bonding, and hence, no relative displacement at the interface was allowed.
Boundary conditions were applied to the RVE for generating strain and stress in the FE model. Nodes in one of two opposing RVE faces were restrained from out-of-plane displacement and rotation, while nodes remained free to move with the plane. This constraint was accomplished defining a cylindrical coordinate system at each respective node with the axial coordinate being aligned in the out-of-plane direction ζ and constraining the appropriate coordinate directions. Nodes in the opposing RVE face were subjected to a uniform displacement of 1 μm perpendicular opposite to the constrained RVE face. As such, a global normal strain εζ of 2500 με was applied to the RVE with a size of 400 μm. Using Equation (5) the modulus of elasticity was computed for each node i located in the RVE face subjected to the prescribed displacement. Boundary conditions were applied to the RVE for generating strain and stress in the FE model. Nodes in one of two opposing RVE faces were restrained from out-of-plane displacement and rotation, while nodes remained free to move with the plane. This constraint was accomplished defining a cylindrical coordinate system at each respective node with the axial coordinate being aligned in the out-of-plane direction ζ and constraining the appropriate coordinate directions. Nodes in the opposing RVE face were subjected to a uniform displacement of 1 µm perpendicular opposite to the constrained RVE face. As such, a global normal strain ε ζ of 2500 µε was applied to the RVE with a size of 400 µm.
Using Equation (5) the modulus of elasticity was computed for each node i located in the RVE face subjected to the prescribed displacement.
where σ ζi and E i are the normal stress extracted from the FE results and modulus of elasticity on the ith node, respectively. As an example, Figure 7 displays normal stress on the model for a FVF of 0.45.
Invoking the Hencky strain method as given by Equation (6) for all nodes subjected to the prescribed displacement yields nodal strains being equivalent to the global strain.
where l i is the deformed RVE length at a specific node position and L the initial RVE length. Finally, the effective modulus of elasticity E eff for the polymer composite embodied by the RVE was calculated using Equation (7).
where n is the total number of nodes.  In order to facilitate evaluating the effect of filler loading on Poisson's ratio, boundary conditions were modified as compared to the modulus of elasticity analysis, that is, Cartesian boundary conditions were applied by restraining the out-of-plane displacement of nodes located in one of two opposing RVE faces, with nodes being free to move in the mutually perpendicular directions. To restrain rigid body motion of the RVE, one corner node of the restrained RVE face was held fixed in all three Cartesian directions. As in the previous case, nodes associated with the RVE face opposite to the constrained plane were again subjected to a 1 μm displacement normal to the plane.
Using Equation (8) In order to facilitate evaluating the effect of filler loading on Poisson's ratio, boundary conditions were modified as compared to the modulus of elasticity analysis, that is, Cartesian boundary conditions were applied by restraining the out-of-plane displacement of nodes located in one of two opposing RVE faces, with nodes being free to move in the mutually perpendicular directions. To restrain rigid body motion of the RVE, one corner node of the restrained RVE face was held fixed in all three Cartesian directions. As in the previous case, nodes associated with the RVE face opposite to the constrained plane were again subjected to a 1 µm displacement normal to the plane.
Using Equation (8), Poisson's ratio values ν i were calculated for each node in the displaced RVE face. The effective Poisson's ratio ν eff for the RVE was determined by Equation (9).
where ε Ψi and ε ζi represent transverse and normal strain, respectively. The effect of filler loading on the CTE was assessed by applying the same constraints to a single RVE face as for the modulus of elasticity analysis. In addition, a temperature change ∆T of 10 K was applied to the entire RVE. The CTE associated with node i, α ζi , and the effective CTE for the RVE, α eff , were calculated using Equations (10) and (11), respectively. As an example, deformations normal to the RVE plane due to thermal expansion are depicted in Figure 8.
where ∆l ζi is axial displacement extracted from the FE results at the ith node.

Results and Discussion
Prior to studying the effect of filler addition on mechanical properties, test cases were performed to assess the distribution of filler particles within RVE. The particle coordinates for the three Cartesian directions were plotted as depicted in the example shown in Figure 9. It can be observed that particles were randomly distributed and dispersed within RVE. Figure 9 also illustrates that, as mentioned above, the largest particles were added to the RVE first (smallest particle index), followed by particles in decreasing size order. Moreover, three sets of numerical analyses were performed, and mechanical properties were calculated for each of the mutually perpendicular axes defining the RVE (X, Y, Z). The example data in Table 2 indicates that differences between the results for given directions are negligible. Hence, it was concluded that the SFEA framework, and especially the RNG module's random number generator, performed satisfactorily for creating true randomness in the analysis.

Results and Discussion
Prior to studying the effect of filler addition on mechanical properties, test cases were performed to assess the distribution of filler particles within RVE. The particle coordinates for the three Cartesian directions were plotted as depicted in the example shown in Figure 9. It can be observed that particles were randomly distributed and dispersed within RVE. Figure 9 also illustrates that, as mentioned above, the largest particles were added to the RVE first (smallest particle index), followed by particles in decreasing size order. Moreover, three sets of numerical analyses were performed, and mechanical properties were calculated for each of the mutually perpendicular axes defining the RVE (X, Y, Z). The example data in Table 2 indicates that differences between the results for given directions are negligible. Hence, it was concluded that the SFEA framework, and especially the RNG module's random number generator, performed satisfactorily for creating true randomness in the analysis.

X-coordinate
Y-coordinate Z-coordinate  After confirming its proper function, the SFEA framework was employed to predict the effective modulus of elasticity, Poisson's ratio, and CTE for the glass bead/epoxy composite for FVF of 0.05, 0.10, 0.15, 0.20, 0.275, 0.35, and 0.45. As mentioned above, the effective properties for each volume fraction were computed via the MCS module and then transferred to the database for statistical analysis. The unbiased standard deviation, variance, and mean were calculated for the stored effective properties data. The MCS module iterates the FEA model until the acceptance criteria defined in the MCS module were satisfied, i.e., the standard deviation reaches an explicit value. Then, the final effective properties for a given FVF were computed as the average from all iterations.
After completing a full Monte Carlo simulation for a specific FVF, a study was performed to assess the quality of the data that was acquired. It was hypothesized that the SFEA framework creates true randomness with data falling under a normal distribution with close congruence between the mean and median values. Ideally, the determination of properties for each FVF would be done with an infinite number of iterations, and hence, predictions were considered continuous random variables for statistical analysis purposes. A probability density function (PDF) was thus computed for each dataset related to a given FVF. The normalized PDF for modulus of elasticity, Poisson's ratio, and CTE data are depicted in Figures 10a, 11a and 12a. The data presented in Figures 10-12 suggests that the effective properties are normally distributed. To test this assumption, statistical data analyses were performed based on mean, median, skewness, and kurtosis values, with skewness and kurtosis providing a measure for asymmetry and peakedness of the distributions, respectively. In other words, these measures were employed to assess how well the data adhere to a normal distribution (e.g., skewness and kurtosis of zero indicate a perfect normal distribution). A comprehensive study was performed by West et al. [46] regarding normal distribution quality indicators. Based on their study, skewness greater than two and kurtosis greater than seven were considered indicators for a significant deviation from normal distribution. In addition, the number of iterations was emphasized as a parameter that had the greatest influence on skewness and kurtosis. Tables 3-5 are showing correspondingly the statistical analyses performed for modeling data relating to the modulus of elasticity, Poisson's ratio, and CTE. The statistical analysis results show that mean and median values are closely congruent, and skewness and kurtosis values are close to zero. Therefore, all datasets were considered to have passed normality test requirements, which is also a confirmation that sufficient numbers of iterations were performed.       As mentioned above, calculated final effective properties can be considered continuous random variables, and hence, the probability of occurrence of a specific value can be calculated within an interval as expressed by Equation (12). Data from the model simulations can thus be presented in the form of cumulative distribution functions (CDFs). CDFs for the effective properties, i.e., the modulus of elasticity, Poisson's ratio, and CTE, are shown in Figures 10b, 11b and 12b, respectively. (12) where P is the probability of an effective property occurring within an interval a and b; X and f (χ) are a continuous random variable and the PDF, respectively. In a final step of the modeling data analysis, results for the effective modulus of elasticity and Poisson's ratio were compared with experimental results published in the technical literature [41,42], as shown in Figures 13 and 14, respectively. For the modulus of elasticity, the experiments indicated a non-linear increase with increasing FVF. It can be observed that the numerical data are in good agreement with the experiments. In fact, the predictions were within the experimental data scatter. Figure 13 also includes predictions from analytical models, namely, the Mori-Tanaka and Hashin-Shtrikman approaches. It can be seen that both analytical methods delivered similar results yet underpredict the values from experiments and numerical modeling, especially for higher FVF.
Referring to Figure 14, experimental data indicated a reduction in Poisson's ratio with rising FVF. The numerical predictions followed the experimental values rather well up to an FVF of 0.15, at which point a shift in the experimental data seemed to have occurred compared to its initial trend. The original publication [42] from which the experimental data were sourced from did not address nor further investigate this behavior. It is herein speculated that at higher FVF, the sample morphology may have deviated from the assumptions made in the present study, which is a random filler distribution and dispersion, leading to a reduced ability to lower the Poisson's ratio.
FVF. The numerical predictions followed the experimental values rather well up to an FVF of 0.15, at which point a shift in the experimental data seemed to have occurred compared to its initial trend. The original publication [42] from which the experimental data were sourced from did not address nor further investigate this behavior. It is herein speculated that at higher FVF, the sample morphology may have deviated from the assumptions made in the present study, which is a random filler distribution and dispersion, leading to a reduced ability to lower the Poisson's ratio.  [41,42]. Error bars express experimental data scatter. Figure 14. Comparison of modeling results for the Poisson's ratio with experimental data [42]. Error bars express experimental data scatter.  [41,42]. Error bars express experimental data scatter.
FVF. The numerical predictions followed the experimental values rather well up to an FVF of 0.15, at which point a shift in the experimental data seemed to have occurred compared to its initial trend. The original publication [42] from which the experimental data were sourced from did not address nor further investigate this behavior. It is herein speculated that at higher FVF, the sample morphology may have deviated from the assumptions made in the present study, which is a random filler distribution and dispersion, leading to a reduced ability to lower the Poisson's ratio.  [41,42]. Error bars express experimental data scatter. Figure 14. Comparison of modeling results for the Poisson's ratio with experimental data [42]. Error bars express experimental data scatter. Figure 14. Comparison of modeling results for the Poisson's ratio with experimental data [42]. Error bars express experimental data scatter.
Numerical predictions for effective CTE are depicted in Figure 15. No relevant experimental data could be found for comparison in the technical literature. Alternatively, a basic Voigt model prediction was made as shown in the figure for comparison with the numerical modeling results. Both the numerical and analytical predictions describe a significant reduction of CTE with increasing FVF. It can be observed that the graph associated with the numerical model was slightly non-linear, predicting a 6% lower CTE for the highest filler loading compared to the Voigt model. Interestingly, the rather basic Voigt model appeared to yield acceptable predictions in this context. However, prior to generalizing such a statement, it would be prudent to study particle size effects and compare data with experimental results.
Both the numerical and analytical predictions describe a significant reduction of CTE with increasing FVF. It can be observed that the graph associated with the numerical model was slightly non-linear, predicting a 6% lower CTE for the highest filler loading compared to the Voigt model. Interestingly, the rather basic Voigt model appeared to yield acceptable predictions in this context. However, prior to generalizing such a statement, it would be prudent to study particle size effects and compare data with experimental results. Similar to other research work, the present study confirms that filler particle and polymer matrix properties, as well as filler loading, significantly affect the effective composite properties. It has been shown that particle size is another important parameter since polymer composites created with smaller particles tend to increase the effective modulus of elasticity to a greater extent compared to composites containing larger particles at the same FVF [47]. Most analytical approaches are unable to capture the influence of filler geometry and, therefore, fail to predict properties accurately [41,48]. Alternatively, empirical models can be fitted to predict the properties of specific composite configurations [49]. While such an approach may increase the accuracy for material property predictions, it has drawbacks in terms of requiring a range of empirical models with limited applicability. The SFEA framework that was applied in the present study alleviates many of the shortcomings associated with analytical and empirical techniques. The SFEA framework, while being computationaly more intensive than analytical and empirical methods, was shown to be versatile in terms of filler and matrix material characteristics and geometries. Moreover, the developed numerical technique is capable of uniting the prediction of mechanical as well as other physical properties [40] in one model, further adding to the versatility of the presented approach. The SFEA framework is thus seen as an attractive tool for material designers. Its application in conjunction with experimental validation may thus serve to accelerate material development to meet application requirements.

Conclusions
A stochastic finite element analysis framework was developed to enable the property prediction of particulate filler modified polymer composites. The employed modeling approach is based on Monte Carlo principles and allows for capturing the effects of filler size distribution, filler shapes, and orientations, which are important parameters for accurately predicting composite properties. Similar to other research work, the present study confirms that filler particle and polymer matrix properties, as well as filler loading, significantly affect the effective composite properties. It has been shown that particle size is another important parameter since polymer composites created with smaller particles tend to increase the effective modulus of elasticity to a greater extent compared to composites containing larger particles at the same FVF [47]. Most analytical approaches are unable to capture the influence of filler geometry and, therefore, fail to predict properties accurately [41,48]. Alternatively, empirical models can be fitted to predict the properties of specific composite configurations [49]. While such an approach may increase the accuracy for material property predictions, it has drawbacks in terms of requiring a range of empirical models with limited applicability. The SFEA framework that was applied in the present study alleviates many of the shortcomings associated with analytical and empirical techniques. The SFEA framework, while being computationaly more intensive than analytical and empirical methods, was shown to be versatile in terms of filler and matrix material characteristics and geometries. Moreover, the developed numerical technique is capable of uniting the prediction of mechanical as well as other physical properties [40] in one model, further adding to the versatility of the presented approach. The SFEA framework is thus seen as an attractive tool for material designers. Its application in conjunction with experimental validation may thus serve to accelerate material development to meet application requirements.

Conclusions
A stochastic finite element analysis framework was developed to enable the property prediction of particulate filler modified polymer composites. The employed modeling approach is based on Monte Carlo principles and allows for capturing the effects of filler size distribution, filler shapes, and orientations, which are important parameters for accurately predicting composite properties. While the developed analysis framework is capable of predicting a variety of mechanical and physical properties, such as thermal conductivity, the present study focused on predicting the effective modulus of elasticity, Poisson's ratio, and coefficient of thermal expansion for the case of randomly distributed and dispersed spherical glass particles embedded in an epoxy polymer matrix.
The numerical results obtained from the analysis framework were compared with experimental values and analytical approaches. The numerical data were found to agree well with non-linear trends observed in the experiments, especially for the elastic modulus, in which case predictions fell within the experimental data scatter. Numerical predictions deviated from experimental Poisson's ratio data for filler volume fractions exceeding 0.15. This may be a result of the morphology changes in test specimens at higher filler loadings, e.g., an increasing extent of filler agglomerations, which would deviate from assumptions made for the numerical model. Nevertheless, in general, the data from the numerical model were found to predict the experiments rather well, while the employed analytical methods were less successful in predicting the test data accurately.
This study shows that numerical techniques, including the current stochastic finite element analysis framework, are attractive and effective for modeling particulate filler modified polymer composites, outperforming analytical methods in terms of versatility and purely experimental campaigns with regard to time and cost. Given that the developed analysis framework is capable of treating binary, ternary, and higher order polymer composites with randomly distributed filler particles, it is seen as a suitable tool for material designers to achieve greater accuracy predictions and use less costly and time-consuming experimentation in order to expedite the development of advanced filler modified polymer composites.
Author Contributions: H.A.M. was chiefly responsible for the conceptualization, methodology, and software development. H.A.M. also performed the formal analysis and validation of results. Following the original draft preparation by H.A.M., further writing of the manuscript, data visualization, review, and editing was supported by P.M., who was also responsible for supervision, project administration, and funding acquisition.