Microstructure-Free Finite Element Modeling for Elasticity Characterization and Design of Fine-Particulate Composites

The microstructure-based finite element modeling (MB-FEM) of material representative volume element (RVE) is a widely used tool in the characterization and design of various composites. However, the MB-FEM has a number of deficiencies, e.g., time-consuming in the generation of a workable geometric model, challenge in achieving high volume-fractions of inclusions, and poor quality of finite element mesh. In this paper, we first demonstrate that for particulate composites the particle inclusions have homogeneous distribution and random orientation, and if the ratio of particle characteristic length to RVE size is adequately small, elastic properties characterized from the RVE are independent of particle shape and size. Based on this fact, we propose a microstructure-free finite element modeling (MF-FEM) approach to eliminate the deficiencies of the MB-FEM. The MF-FEM first generates a uniform mesh of brick elements for the RVE, and then a number of the elements, with their total volume determined by the desired volume fraction of inclusions, is randomly selected and assigned with the material properties of the inclusions; the rest of the elements are set to have the material properties of the matrix. Numerical comparison showed that the MF-FEM has a similar accuracy as the MB-FEM in the predicted properties. The MF-FEM was validated against experimental data reported in the literature and compared with the widely used micromechanical models. The results show that for a composite with small contrast of phase properties, the MF-FEM has excellent agreement with both the experimental data and the micromechanical models. However, for a composite that has large contrast of phase properties and high volume-fraction of inclusions, there exist significant differences between the MF-FEM and the micromechanical models. The proposed MF-FEM may become a more effective tool than the MB-FEM for material engineers to design novel composites.


Introduction
Particulate composites are widely used in industrial products and engineering structures due to their merits, such as ease-of-manufacturing and great design flexibility. At the length scale of material representative volume element (RVE), the behavior of a particulate composite can be usually considered as homogeneous and isotropic if (1) the geometric aspect ratio of inclusions is small, (2) both the distribution and orientation of inclusions are statistically homogeneous, and (3) the characteristic dimension of inclusions is adequately small compared with the size of RVE. An important task in the design of composites is to predict their elastic properties based on phase properties and volume fractions. A number of methods is available for the prediction, and they are often classified into experimental, analytical, and numerical categories. Experimental methods are direct and reliable, but they are also expensive and time consuming. Analytical solutions developed from micromechanics models are convenient and efficient; however, they have various limitations in application because they are based on special assumptions regarding composite microstructure. Numerical methods, mainly represented by finite element modeling of RVE, 2 of 14 have a number of advantages over the other two categories of methods. They are more efficient than experimental methods, and they do not make any special assumption about composite microstructure. The microstructure-based finite element modeling (MB-FEM) of material RVE has played an important role for material engineers to understand the micromechanics of composite cracking and debonding, which often initiate and develop in the phase materials or at their interfaces [1][2][3]. An accurate representation of the composite microstructure is critical in the finite element analysis of microscopic damage.
However, for the design of a novel composite, the main interest is in the elastic properties at the macroscopic level, rather than the interactions of phase materials at the microscopic level. The commonly considered design variables are elastic properties and volume fractions of phase materials. The design process is often iterative and involves a large number of finite element analyses of intermediate designs; the conventional MB-FEM has a number of deficiencies for this purpose. First, the creation of a workable geometric model for the RVE microstructure is often tedious and time-consuming. A number of algorithms has been developed to generate RVE of particulate composites, e.g., [4][5][6][7][8][9][10][11][12][13] among others. These algorithms either have difficulty to accommodate a high volumefraction of inclusions [7][8][9]13], or have to use a time-consuming iterative process to achieve a high volume-fraction with the sacrifice of microstructure randomness [10][11][12]. It should be mentioned that the majority of the algorithms are developed for simple inclusion shapes such as sphere, ellipsoid, and cylinder. Inclusion shapes in actual composite materials are often irregular, which would further complicate the creation of valid RVE geometric models. Geometric models of such composites usually contain a large number of degenerated and small geometric entities. Therefore, even with the success of creating an RVE geometric model, there are still challenges in the generation of valid and high-quality finite element meshes, because the degenerated and small entities in the geometric models either prevent the generation of finite element meshes or result in poor-quality elements. Iterative mesh adaptation algorithms usually have to be used to improve mesh quality, but a satisfactory quality is still not guaranteed.
To eliminate the limitations of the conventional MB-FEM for the characterization and design of fine-particulate composites, we first demonstrate that if the size ratio of inclusion to RVE is adequately small, elastic properties of RVE computed by finite element modeling are independent of inclusion shape and size. Based on this fact, we propose microstructure-free finite element modeling (MF-FEM). We validate the MF-FEM against experimental data reported in the literature and compare the MF-FEM with the widely used micromechanical models.

Effect of Inclusion Shape and Size on Elastic Properties of Particulate Composites
To study the effect of inclusion shape and size on the elastic properties of particulate composites, the following assumptions are made: (1) the materials of the matrix and the inclusions are homogeneous, isotropic, and perfectly bond to each other; (2) the distribution and orientation of inclusions of small aspect ratio are statistically homogeneous; (3) the applied loading only introduces elastic deformation. Under the above assumptions, we hypothesize that for a fixed volume-fraction of inclusions, if the ratio of inclusion characteristic length to RVE size is adequately small, composite elastic properties such as Young's modulus, Poisson's ratio, and shear modulus are independent of inclusion shape and size, and also independent of loading orientation. To verify the above hypothesis, a series of particulate-composite RVE models, which are different from each other only by their inclusion shape and size, are constructed and analyzed. The models are described as follows: • All RVE models are in a cubic shape, and the sides have a length of 100 units. Since inclusion-to-RVE size ratio is of interest, the unit can be in any length from nanometer to meter.

•
The composite is a particulate-filled glassy polymer [14]. The matrix material has Young's modulus E m = 2.68 GPa and Poisson's ratio V m = 0.394, and the inclusion material has Young's modulus E i = 70.0 GPa and Poisson's ratio V i = 0.23.

•
The volume fraction of inclusions in all models is fixed at 30%.

•
The inclusions in a model have the same shape, i.e., either spheroid, almond-shaped, or pill-shaped.

•
The inclusions in a model have the same size. The inclusion-to-RVE size ratio is gradually reduced from 0.8 to 0.02 in the series of models, by decreasing inclusion sizes.
A set of such models is shown in Figure 1.
i. 2022, 6, x FOR PEER REVIEW 3 of 13 • All RVE models are in a cubic shape, and the sides have a length of 100 units. Since inclusion-to-RVE size ratio is of interest, the unit can be in any length from nanometer to meter.

•
The composite is a particulate-filled glassy polymer [14]. The matrix material has Young's modulus = 2.68 GPa and Poisson's ratio = 0.394, and the inclusion material has Young's modulus = 70.0 GPa and Poisson's ratio = 0.23.

•
The volume fraction of inclusions in all models is fixed at 30%.

•
The inclusions in a model have the same shape, i.e., either spheroid, almond-shaped, or pill-shaped.

•
The inclusions in a model have the same size. The inclusion-to-RVE size ratio is gradually reduced from 0.8 to 0.02 in the series of models, by decreasing inclusion sizes.
A set of such models is shown in Figure 1.  RVE Young's modulus, Poisson's ratio, and shear modulus are characterized by the conventional MB-FEM. The coordinate system used in the finite element modeling is illustrated in Figure 2. Boundary conditions for characterizing the elastic properties are listed in Table 1. RVE Young's modulus, Poisson's ratio, a conventional MB-FEM. The coordinate system lustrated in Figure 2. Boundary conditions f listed in Table 1.

RVE Surface
Young's Modulus (E i , i = x, y, z) and Poisson's Ratio (ν ij , i, j = x, y, z) Shear Modulus (G ij , i, j = x, y, z) Based on the mean-field homogenization theory [15], RVE Young's modulus (E i ), shear modulus (G ij ), and Poisson's ratio (ν ij ) are determined from RVE average stresses (σ i , τ ij ) and average strains (ε i , γ ij ): The average stresses and strains are calculated as In Equation (4), σ i , ε i , τ ij , and γ ij are, respectively, the normal stress, normal strain, shear stress, and shear strain that are determined by the finite element modeling; V is the volume of the RVE.
All finite element analyses in this study were conducted using commercial software, ANSYS Mechanical APDL (2020 R1). For each model, the mean value (P) and standard deviation (σ) of Young's moduli, shear moduli, and Poisson's ratios characterized from the three loading orientations were calculated, i.e., where P represents one of the three elastic properties. Variations of the mean values and standard deviations with the inclusion-to-RVE size ratio are shown in Figure 3.
The following observations can be made from the results shown in Figure 3: • With the inclusion-to-RVE size ratio approximately smaller than 0.04, the RVE elastic properties are almost not affected by the inclusion shape and size. • Anisotropy in the RVE elastic properties, which is measured by the error bars in Figure 3, is gradually reduced and then disappears with the decreased inclusion-to-RVE size ratio.
It is not surprising that the hypothesis is valid, because it is the fundamental assumption adopted in the continuum mechanics of composite materials, and it has been demonstrated to be true in a large number of previous studies, which were conducted on different inclusion shapes and phase properties, e.g., [16][17][18][19][20][21][22][23] among others. The threshold of inclusion-to-RVE size ratio referred in the above is equivalent to the critical RVE size used in other studies. Trias et al. [21] reported that for a typical unidirectionally carbonfiber-reinforced polymer, to properly represent the random distribution of fibers in the transverse plane, the minimum size of RVE should be 50 times the fiber radius. Harper et al. [16] showed that for a planar model of discontinuous carbon fiber composites, the convergence of results from multiple realizations (i.e., models of the same volume fraction) is achieved at RVE edge lengths that are four times the fiber length, irrespective of fiber volume fraction; however, if only one realization is allowed, the ratio of fiber length to RVE size must be reduced to 0.03. The studies show that the averaging of multiple realizations has the effect of reducing anisotropy and thus underestimating the effect of RVE size or inclusion-to-RVE size ratio. For this reason, averaging of multiple realizations of the same volume fraction is not used in this study, and the results in Figure 3 were produced from one realization. It should be pointed out that the critical size of particulate-composite RVE can be significantly different in different regimes, for example, linear vs. nonlinear properties [22]. where represents one of the three elastic properties. Variations of the mean valu standard deviations with the inclusion-to-RVE size ratio are shown in Figure 3.

Microstructure-Free Finite Element Modeling of Particulate-Composites
The validity of the hypothesis demonstrated in the previous section justifies the idea of microstructure-free finite element modeling (MF-FEM) for particulate composites. Now that RVE properties are independent of inclusion shape and size if the inclusion-to-RVE size ratio is adequately small, theoretically we can consider any inclusion shape and size that are convenient for the finite element modeling. In this section, we propose a microstructure-free finite element modeling procedure for the characterization and design of two-phase particulate composites. The procedure is similar to the conventional MB-FEM. The only difference is in the representation of inclusions. A geometric model of the RVE microstructure is not needed; instead, the inclusions are represented by finite elements. The procedure is described below.

•
A cubic RVE with side length of 100 units is constructed. The unit can be at any length scale from nanometer to meter, depending on the composite material to be studied.

•
The RVE is meshed with brick elements of the same size, which is determined by the critical inclusion-to-RVE size ratio, i.e., 0.04 to 0.02 times the RVE side length, as discussed in the previous section. All the elements are first assigned with the properties of the matrix material. • Then, a number of the elements is randomly selected and re-assigned with the properties of the inclusion material. The number of the selected elements is determined by the desired volume fraction of inclusions and the volume of each element. Samples of such microstructure-free finite element models are shown in Figure 4. • Then, the boundary conditions described in Table 1   The effect of element-to-RVE size ratio on elastic properties predicted by the MF-FEM was studied using the same composite described in the previous section, with the volume fraction of inclusion elements set to 30%. The variations of RVE properties, averaged in the three loading orientations, with element-to-RVE size ratio, are shown in Figure 5. Similar to those computed by the MB-FEM, RVE properties predicted by the MF-FEM also show anisotropy if a large element-to-RVE size ratio is used; the degree of anisotropy is indicated by the error-bar length in the figure. With the element-to-RVE size ratio decreased, the anisotropy disappears, as shown in Figure 5. By comparing Figures 3 and 5, one can observe that the variations of properties with element-to-RVE size ratio shown in Figure 5 are much smaller than those in Figure 3, indicating that even with a relatively larger element-to-RVE size ratio, and thus a smaller number of elements, the MF-FEM is The effect of element-to-RVE size ratio on elastic properties predicted by the MF-FEM was studied using the same composite described in the previous section, with the volume fraction of inclusion elements set to 30%. The variations of RVE properties, averaged in the three loading orientations, with element-to-RVE size ratio, are shown in Figure 5. Similar to those computed by the MB-FEM, RVE properties predicted by the MF-FEM also show anisotropy if a large element-to-RVE size ratio is used; the degree of anisotropy is indicated by the error-bar length in the figure. With the element-to-RVE size ratio decreased, the anisotropy disappears, as shown in Figure 5. By comparing Figures 3 and 5, one can observe that the variations of properties with element-to-RVE size ratio shown in Figure 5 are much smaller than those in Figure 3, indicating that even with a relatively larger element-to-RVE size ratio, and thus a smaller number of elements, the MF-FEM is still able to predict RVE properties with reasonable accuracy. the three loading orientations, with element-to-RVE size ratio, are shown in Figure 5. Similar to those computed by the MB-FEM, RVE properties predicted by the MF-FEM also show anisotropy if a large element-to-RVE size ratio is used; the degree of anisotropy is indicated by the error-bar length in the figure. With the element-to-RVE size ratio decreased, the anisotropy disappears, as shown in Figure 5. By comparing Figures 3 and 5, one can observe that the variations of properties with element-to-RVE size ratio shown in Figure 5 are much smaller than those in Figure 3, indicating that even with a relatively larger element-to-RVE size ratio, and thus a smaller number of elements, the MF-FEM is still able to predict RVE properties with reasonable accuracy. Comparison of RVE properties computed by the MF-FEM and the MB-FEM is displayed in Figure 6. With the inclusion-to-RVE (or element-to-RVE for the MF-FEM) size ratio of 0.02, differences in RVE properties computed by the two types of finite element models are negligible, suggesting that the MF-FEM is able to replace the MB-FEM in the Comparison of RVE properties computed by the MF-FEM and the MB-FEM is displayed in Figure 6. With the inclusion-to-RVE (or element-to-RVE for the MF-FEM) size ratio of 0.02, differences in RVE properties computed by the two types of finite element models are negligible, suggesting that the MF-FEM is able to replace the MB-FEM in the prediction of effective properties of particulate composites.

Validation of MF-FEM against Experimental Data and Comparison with Popular Micromechanics Models
In this section, the MF-FEM was validated against experimental data previously reported in the literature and compared with popular micromechanics models. Two repre-

Validation of MF-FEM against Experimental Data and Comparison with Popular Micromechanics Models
In this section, the MF-FEM was validated against experimental data previously reported in the literature and compared with popular micromechanics models. Two representative cases of experimental data were selected for the validation. One is WCcobalt alloy, where the phase properties are close to each other [24], and the other is glass particulate reinforced polyester, where the phase properties have a large contrast [25]. Using a dynamic resonance method, Doi et al. [24] measured the Young's modulus, shear modulus, bulk modulus, and Poisson's ratio for a series of WC-cobalt alloys, with volume fractions of WC varying from 0.55 to 0.91. Properties of WC and cobalt used in the study are listed in Table 2. Richard [25] conducted a series of tensile tests on specimens of glass microsphere filled polyester composites to measure their Young's modulus and Poisson's ratio; the composites had volume fractions of glass microsphere varying from 0.0 to 0.464. The phase material properties used in Richard's study are listed in Table 3. Bulk modulus and shear modulus of the composites and phase materials were calculated using the elasticity relationships. Although there exists a large number of micromechanics models available for comparison [23,26,27], we selected the models that have explicit analytical solutions, explicitly developed for particulate composites and widely used in the current literature. For the glass-polyester composites, as can be seen from Figure 8a,b, at low volume fractions, the MF-FEM predictions were well congruous with the experimental data and the predictions of GSC and M-T methods. However, a deviation trend occurred among the curves with the increase of glass volume fraction. Unfortunately, experimental data of high volume-fraction of glass were not available for the validation. Experimental data of bulk modulus and Poisson's ratio, as displayed in Figure 8c,d, showed much greater scatterness, and the agreement between the experimental data and the MF-FEM predictions was not as good as in the Young's modulus or shear modulus.
It should be pointed out that the specific loading and constraint conditions used in the experiments are not available in the literature, which may be considerably different from those described in Table 1 and used in the MF-FEM. This may explain the differences between the experimental data and the MF-FEM predictions. From Figures 7 and 8 it can be noticed that compared with the experimental data, micromechanics models such as GSC and M-T underestimate composite Young's modulus, shear modulus, and bulk modulus, which may become even more significant for composites with a large contrast of phase properties and a high volume-fraction of inclusions. On the other hand, the predictions of MF-FEM are generally higher than the experimental data; the possible reason is that porosity may exist in the experimental materials, while the MF-FEM considers the materials as non-porous.
Glass microsphere 68.95 0.21 Although there exists a large number of micromechanics models available for comparison [23,26,27], we selected the models that have explicit analytical solutions, explicitly developed for particulate composites and widely used in the current literature. Based on the above considerations, we chose Voigt-Reuss (V-R) bounds, Hashin-Shtrikman (H-S) bounds, the generalized self-consistent (GSC) scheme, the Mori-Tanaka (M-T) method, and the Voigt-Reuss-Hill (V-R-H) average for comparison. Figures 7 and 8 show the results of validation and comparison. As it can be observed from Figure 7, for the WC-cobalt alloys, all the properties computed by the MF-FEM, except Poisson's ratio, had excellent agreement with the experimental data and the predictions by the GSC and the M-T methods.

Discussion
Although there exist differences between the MF-FEM predictions and the experimental data, composite properties computed by the MF-FEM are in excellent agreement with those produced by the MB-FEM, as shown in Figure 6. The MF-FEM will become a more powerful and more efficient tool than the conventional MB-FEM in the design of particulate composites.
The MF-FEM has a number of advantages over the MB-FEM for the prediction of composite properties. There is no need to create a geometric model for the composite microstructure, which is usually tedious and time-consuming. There is no need to determine the size of RVE. Mesh quality is always guaranteed, and the number of elements does not change with the complexity of the composite microstructure. Once a universal element-to-RVE size ratio is determined, even a convergence study is not required, which is always necessary for the conventional MB-FEM. A high volume-fraction of inclusion can be easily achieved in the MF-FEM.
However, limitations of the MF-FEM are also obvious. The MF-FEM is not able to predict stress distribution in the phase materials, because the microstructure is not represented. If the inclusion-to-RVE size ratio is not adequately small, composite properties predicted by the MF-FEM will have low accuracy.

Conclusions
Based on the study results, it can be reasonably concluded that the MF-FEM is able to replace the conventional MB-FEM for the prediction of composite properties at the macroscopic level, but the MF-FEM is not effective for the study of composite damage at the microscopic level.