International Journal of Composite Materials

p-ISSN: 2166-479X    e-ISSN: 2166-4919

2021;  11(3): 54-61


Received: Oct. 4, 2021; Accepted: Nov. 1, 2021; Published: Nov. 15, 2021


A Numerical Comparison between Geometry and Voxel Based Finite Element Modeling for Elasticity Characterization of Particulate Composites

Yunhua Luo

Department of Mechanical Engineering, University of Manitoba, Winnipeg, R3T 2N2, Canada

Correspondence to: Yunhua Luo, Department of Mechanical Engineering, University of Manitoba, Winnipeg, R3T 2N2, Canada.


Copyright © 2021 The Author(s). Published by Scientific & Academic Publishing.

This work is licensed under the Creative Commons Attribution International License (CC BY).


Elasticity prediction based on phase properties is an important task in the analysis and design of composite materials. Microstructure-based finite element modeling of composite representative volume element (RVE) has a number of advantages over experimental and analytical methods for the task. However, there are a number of challenges in the conventional geometry-based finite element modeling (GB-FEM) of RVE, e.g. time-consuming in the construction of workable geometric models and in the generation of high-quality finite element meshes. To overcome these challenges, we developed a voxel-based finite element modeling (VB-FEM) procedure. In contrast to GB-FEM, VB-FEM first generates a uniform grid mesh and then identifies elements that belong to the inclusions. The rest steps are the same as those in GB-FEM. We compared the performance of GB-FEM and VB-FEM in two representative numerical examples. The results show that upon convergence, elasticity constants characterized by VB-FEM and GB-FEM have excellent agreement. Regarding their computational efficiency, although for the composite having regularly-distributed and large-size inclusions, GB-FEM is significantly more efficient than VB-FEM; However, for the composite having randomly-distributed and small-size inclusions, VB-FEM has similar and even higher efficiency than GB-FEM. Besides overcoming the mentioned challenges, VB-FEM has a number of additional advantages over GB-FEM. We conclude that VB-FEM is an effective tool for elasticity prediction of particulate composites.

Keywords: Particulate composite, Composite elastic property, Representative volume element, Finite element modeling, Geometry-based, Voxel-based

Cite this paper: Yunhua Luo, A Numerical Comparison between Geometry and Voxel Based Finite Element Modeling for Elasticity Characterization of Particulate Composites, International Journal of Composite Materials, Vol. 11 No. 3, 2021, pp. 54-61. doi: 10.5923/j.cmaterials.20211103.02.

1. Introduction

Particulate composite materials are widely used in industrial products and engineering structures due to their merits such as easy-to-manufacture, customized mechanical properties and great flexibility of design. Prediction of composite elasticity constants based on phase properties is a key task in the design of various composite materials. There are primarily three categories of methods available for the prediction of composite properties, i.e. mechanical testing [1-3], analytical micromechanics models [4,5] and finite element modeling of representative volume element (RVE) [6-9]. Mechanical testing is direct and reliable, but it is also time-consuming and expensive. Micromechanical models are the most efficient and convenient ones, because they offer analytical expressions. Nevertheless, they are usually based on special assumptions, for example dilute dispersion. Therefore, they have limited accuracy for composites that do not strictly satisfy the assumptions, e.g. composites that have large contrast of phase properties and high volume-fraction of inclusions [5]. Finite element modeling of RVE has the merits of both mechanical testing and micromechanical models, but avoids their drawbacks. Therefore, finite element modeling of RVE is widely applied in the study and design of composite materials [10-14].
There are basically two types of finite element modeling approaches available for elasticity prediction of composite materials, one is geometry based [10,15], the other is voxel based [16,17]. Take particulate composite as an example, in geometry-based finite element modeling (GB-FEM), a geometric model of the microstructure of RVE, which includes matrix and inclusions, is first constructed; Material properties are assigned to the matrix and the inclusions; A finite element mesh is then generated from the geometric model; Appropriate loading and constraint conditions are applied for the characterization of various elastic properties. Although theoretically geometry-based finite element modeling is an ideal tool for composite property prediction, there are a number of practical challenges. The first one is in the construction of a workable RVE geometric model [10,18]. Inclusions in actual composites are random distributed and have different geometric shapes and sizes. Inclusions (e.g. two spheres) that are infinitely close to each other may result in severely distorted and even degenerated geometric entities. This kind of geometric model is difficult to create even with relaxed geometric tolerances. Even with a valid RVE geometric model, there are still challenges in the generation of high-quality finite element meshes [19,20]. The existence of distorted and degenerated geometric entities may first fail the generation of a finite element mesh and then may result in severely distorted finite elements. Mesh quality significantly affects the accuracy of finite element solutions. To improve mesh quality, adaptive meshing has to be used [20], which is time-consuming but still cannot guarantee the quality.
Voxel-based finite element models (VB-FEM) constructed from quantitative computed tomography (CT) or micro-CT are often used in the study of bone mechanical properties [17,21-23]. The main difference between VB-FEM and GB-FEM is in the generation of finite element mesh. In contrast to GB-FEM, where a finite element mesh is generated from a geometric model, VB-FEM directly converts voxels in digital images into brick elements. Another advantage of VB-FEM is that material properties are reflected in voxel intensity or Hounsfield Unit (HU), which provides realistic information on the distribution of inhomogeneous bone density. There are two major concerns about VB-FEM. One is its computational efficiency. Since digital images usually consist of a large number of voxels, the resulting finite element meshes usually consist of a huge number of elements. The other concern is the non-smoothness of model boundaries, which may introduce stress concentration.
Despite of the above concerns, VB-FEM avoids the difficulties of GB-FEM in the creation of geometric model and in the generation of high-quality mesh. Voxel-based finite element meshes always have good quality. For the above reasons, we are motivated to compare the accuracy and efficiency of GB-FEM and VB-FEM in elasticity prediction of particulate composites, to study the possibility of using VB-FEM in elasticity prediction of particulate composites to avoid the difficulties of GB-FEM.

2. Materials and Methods

To compare their performance, GB-FEM and VB-FEM are constructed from the same composite RVE as illustrated in Figure 1. In GB-FEM, the matrix and the inclusions are meshed separately. A number of issues are introduced by this strategy, which has been discussed before. In contrast, for VB-FEM a uniform grid mesh is first generated for the whole domain occupied by the RVE. Then, the information of inclusion location and size in the RVE are used to identify the elements that belong to the matrix and the inclusions respectively. The resulting finite element mesh with distinguished Table 1 matrix and inclusions looks like a binary digital image, therefore, the terminology of voxel is used. As it will be shown later, the resolution of VB-FEM mesh, defined as the ratio between element size and the RVE size, affects the accuracy of VB-FEM. The rest steps of finite element analysis, including assignment of material properties and application of loading and constraint conditions, are the same as in GB-FEM. The objective of this comparison study is to find out whether or not GB-FEM and VB-FEM have good agreement upon their convergence, and what are their relative computational efficiency.
Table 1. RVE boundary conditions for characterization of composite elastic properties
Figure 1. Construction of GB-FEM and VB-FEM from the same composite RVE
RVE elastic properties, including Young’s modulus, Poisson’s ratio and shear modulus, are predicted by GB-FEM and VB-FEM using the same loading and constraint conditions. A cubic RVE with the setup of a coordinate system is shown in Figure 2. The side length of the RVE is L = 100 units, the unit can be at any length-scale from nanometer to meter, because in principle both GB-FEM and VB-FEM are applicable to problems of any length-scale.
Figure 2. Cubic RVE and coordinate system
The loading and constraint conditions listed in Table 1 are applied to the RVE for the prediction of Young’s modulus , Poisson’s ratio and shear modulus . The applied constraint conditions are intended to mimic the behavior of the RVE in a larger material body and to avoid ‘artificial’ stress concentration. Take Young’s modulus as an example, the RVE surfaces at x = 0, y = 0 and z =0 are constrained in x, y and z direction, respectively; a uniform displacement is introduced on surface x = L; the surfaces y = L and z = L are forced to have uniform displacement from Poisson’s effect.
After the solution of the finite element equations, RVE Young’s modulus, Poisson’s ratio and shear modulus are determined by the homogenization theory based on volume-averaged stresses and strains [24-26], i.e.
In the above equations, , and are respectively the normal stress, normal strain, shear stress and shear strain determined by the finite element analysis; V is the volume of RVE.
Agreement between GB-FEM and VB-FEM is measured by the following normalized difference,
Where P is a property of the RVE, for example, Young’s modulus, Poisson’s ratio, and shear modulus; is the elastic property predicted by GB-FEM upon convergence; is the property obtained by VB-FEM with a certain mesh resolution. A smaller δ value indicates a better agreement between GB-FEM and VB-FEM; If δ = 0, there would be no difference.
To compare computational efficiency of GB-FEM and VB-FEM, normalized number of nodes (η) is used, which is defined as
Where and are respectively the number of finite element nodes in VB-FEM and GB-FEM. If η > 0, more nodes are used in VB-FEM than in GB-FEM; otherwise less nodes are used.

3. Numerical Examples and Results

A polymer filled with small glass spheres [3] is selected for this comparison study. The polymer (matrix) has Young’s modulus = 2.68 GPa and Poisson’s ratio = 0.394; The glass spheres (inclusions) have Young’s modulus = 70.0 Gpa and Poisson’s ratio = 0.23. Two types of RVE of the composite, as shown in Figure 3 and Figure 4, are constructed to compare GB-FEM and VB-FEM. RVE 1 in Figure 3 has large size and regularly distributed spheres, while RVE 2 in Figure 4 has small size and randomly distributed spheres.
Figure 3. RVE 1: regularly-distributed and large-size spheres. (a) composite; (b) matrix; (c) inclusions
Figure 4. RVE 2: randomly-distributed and large-size spheres. (a) composite; (b) matrix; (c) inclusions
All finite element results are produced using ANSYS Mechanical APDL (2020 R1, ANSYS, Inc.). In GB-FEM, 6-node tetrahedral element is used with mesh adaptation; In VB-FEM, 8-node brick element is used.
For RVE 1, sample GB-FEM and VB-FEM meshes are displayed respectively in Figure 5 and Figure 6. Variations of Young’s moduli predicted by GB-FEM and VB-FEM with the number of element nodes are shown in Figure 7.
Figure 5. GB-FEM of RVE 1. (a) whole model; (b) matrix; (c) inclusions
Figure 6. VB-FEM of RVE 1. (a) whole model; (b) matrix; (c) inclusions
Figure 7. Variation of RVE Young’s modulus with number of nodes (RVE 1)
Normalized differences between GB-FEM and VB-FEM are provided in Table 2, where the normalized differences (δ) are calculated based on Equation (5) and the normalized number of nodes is determined by Equation (6). It should be men-tioned that δ(V) in Table 2 (as well in Table 3) is calculated from the volumes of elements representing the spherical inclusions, there is a difference between the sum of the element volumes and that of the spheres.
Table 2. Normalized differences (%) between VB-FEM and converged GB-FEM (RVE 1)
Table 3. Normalized differences (%) between VB-FEM and converged GB-FEM (RVE 2)
For RVE 2, sample meshes are displayed in Figure 8 and Figure 9. Variations of Young’s moduli predicted by GB-FEM and VB-FEM with the number of element nodes are shown in Figure 10, the Young’s modulus in the figure is the average of Young’s moduli in the x, y and z axial direction shown in Figure 2. The normalized differences between GB-FEM and VB-FEM are provided in Table 3.
Figure 8. GB-FEM of RVE 2. (a) whole model; (b) matrix; (c) inclusions
Figure 9. VB-FEM of RVE 2. (a) whole model; (b) matrix; (c) inclusions
Figure 10. Variation of RVE Young’s modulus with number of nodes (RVE 2)
The following observations can be made from the results:
1) With mesh resolution refined in VB-FEM, the normalized difference in volume fraction of inclusion, i.e. δ(V) in Table 2 and Table 3, is gradually reduced, indicating that VB-FEM and GB-FEM have closer volume fraction of inclusions, but it is difficult for them to have exactly the same.
2) VB-FEM and GB-FEM converge to very similar final solutions, which can be observed from the convergence curves in Figure 7 and Figure 10, and also from the normalized differences in Table 2 and Table 3. If the differences in the volume fraction of inclusions are taken in account, the solutions of VB-FEM and GB-FEM would be even closer.
3) Even with a coarse mesh resolution, VB-FEM solutions have good agreement with the converged GB-FEM. For example, with a 1:20 mesh resolution, the maximum normalized difference is only 1.10% for RVE 1 (Table 2) and 2.36% for RVE 2 (Table 3). Elimination of the differences in volume fraction of inclusions would further improve the agreement.
4) Regarding computational efficiency, GB-FEM has much greater advantage over VB-FEM for RVE 1, which can be seen from the convergence in Figure 7 and the normalized number of element nodes (η) in Table 2. However, for RVE 2, VB-FEM and GB-FEM have very similar efficiency, as shown by the convergence curve in Figure 10 and the η values in Table 3.

4. Discussion and Conclusions

The results of this comparison study show that for the elasticity prediction of particulate composites, GB-FEM and VB-FEM have similar accuracy upon their convergence. Regarding their computational efficiency, for a composite with large size and regularly distributed inclusions (e.g. RVE 1), GB-FEM has obvious advantage over VB-FEM; but for a composite with small size and randomly distributed inclusions (e.g. RVE 2), which is more often encountered in engineering, VB-FEM and GB-FEM have similar computational efficiency. Nevertheless, VB-FEM avoids all the difficulties of GB-FEM that have been discussed at the beginning of this paper. In addition, VB-FEM has a number of other advantages over GB-FEM. For the creation of VB-FEM, the only requirement is nonoverlap among inclusions, which is much easier to satisfy than the creation of a workable geometric model in GB-FEM. There is no restriction on how close two inclusions can be in VB-FEM; therefore, a high volume-fraction of inclusion can be achieved more easily. The accuracy and computational efficiency of VB-FEM are not affected by the number and the size of inclusions in the composite, as can be seen from the results of RVE 1 and RVE 2 obtained by VB-FEM. However, GB-FEM has to greatly increase the number of elements, because a large number of small-size elements are required to realistically represent the small inclusions and to fill in the very narrow gaps between infinitely close inclusions. If the volume fraction of inclusions is high, the situation will become even worse for GB- FEM [27].
Although GB-FEM and VB-FEM have similar accuracy in the prediction of composite elastic properties, there exist significant differences in the stresses and strains determined by the two models. Figure 11 (a) and (b) show the patterns of von-Mises (v-M) stress distribution produced by GB-FEM and VB-FEM for RVE 1 respectively. Quantitative differences can be observed from the results listed in Table 4.
Table 4. Maximum stresses and strains in composite matrix and inclusion (RVE 1)
Figure 11. Von-Mises stress distribution in RVE 1. (a) GB-FEM; (b) VB- FEM
The results in Table 4 show that, generally, with the mesh resolution refined, the maximum stresses and maximum strains predicted by VB-FEM gradually approach those by GB-FEM. However, even with a fine resolution such as 1:100, VB-FEM still overestimates the stresses and strains compared with GB-FEM. The main cause is that the interface between inclusion and matrix is not smooth in VB-FEM, which may cause stress concentration at the interface. It is well known that stresses and strains predicted by GB-FEM are also sensitive to mesh quality, especially in composites with irregularly distributed inclusions. As discussed before, distorted and degenerated geometric entities in GB-FEM produce distorted elements, which introduce spuriously large stresses and strains.
Based on this study, it can be concluded that VB-FEM can be used in replace of GB-FEM to predict elastic properties of particulate composites, but probably not suitable for predicting stresses and strains within the composites.


The reported research has been supported by Natural Sciences and Engineering Research Council (NSERC) of Canada in grant RGPIN-2019-05372, which is gratefully acknowledged.


[1]  H. Doi, Y. Fujiwara, K. Miyake, and Y. Oosawa. A systematic investigation of elastic moduli of WC-Co alloys. Metallurgical and Materials Transactions B, 1:1417 – 1425, 1970.
[2]  T.G. Richard. The mechanical behavior of a solid microsphere filled composite. J. Composite Materials, 9:108 – 113, 1975.
[3]  J.C. Smith. Experimental values for the elastic constants of a particulate-filled glassy polymer. J. Res. Natl. Bur. Stand. A: Phys. Chem., 80:45 – 49, 1976.
[4]  B. Raju, S.R. Hiremath, and D. Roy Mahapatra. A review of micromechanics based models for effective elastic properties of reinforced polymer matrix composites. Composite Structures, 204:607 – 619, 2018.
[5]  R.M. Christensen. A critical evaluation for a class of micromechanics models. J. Mech, Phys. Solids, 38:379 – 404, 1990.
[6]  N. Chawla, B.V. Patel, M. Koopman, K.K. Chawla, R. Saha, B.R. Patterson, E.R. Fuller, and S.A. Langer. Microstructurebased simulation of thermomechanical behavior of composite materials by object-oriented finite element analysis. Materials Characterization, 49:395 – 407, 2002.
[7]  H.K. Park, J. Jung, and H.S. Kim. Three-dimensional microstructure modeling of particulate composites using statistical synthetic structure and its thermo-mechanical finite element analysis. Computational Materials Science, 126:265 – 271, 2017.
[8]  B. Máša, L. Náhlík, and P. Hutar. Particulate composite materials: Numerical modeling of a cross-linked polymer reinforced with alumina-based particles. Mechanics of Composite Materials, 49:421 – 428, 2013.
[9]  I.V. Singh, A.S. Shedbale, and B.K. Mishra. Material property evaluation of particle reinforced composites using finite element approach. Journal of composite materials, 50:2757 – 2771, 2016.
[10]  S. Bargmann, B. Klusemann, J. Markmannc, J.E. Schnabel, K. Schneider, C. Soyarslan, and J. Wilmers. Generation of 3D representative volume elements for heterogeneous materials: A review. Progress in Materials Science, 96: 322– 384, 2018.
[11]  E.P. David Müzel, S. and; Bonhin, N.M. Guimarães, and E.S. Guidi. Application of the finite element method in the analysis of composite materials: A review. Polymers, 12:818, 2020.
[12]  P. Kenesei, A. Borbély, and H. Biermann. Microstructure based three-dimensional finite element modeling of particulate reinforced metal–matrix composites. Materials Science and Engineering: A, 387 - 389:852 – 856, 2004.
[13]  V. Kukshal, S. Gangwar, and A. Patnaik. Experimental and fi- nite element analysis of mechanical and fracture behavior of SiC particulate filled A356 alloy composites: Part I. Journal of Materials Design and Applications, 229:91 – 105, 2013.
[14]  J. Zhang, Q. Ouyang, Q. Guo, and et al. 3D Microstructure- based finite element modeling of deformation and fracture of SiCp/Al composites. Composites Science and Technology, 123:1 – 9, 2016.
[15]  G. Jagadeesh and S. Gangi Setti. A review on micromechanical methods for evaluation of mechanical behavior of particulate reinforced metal matrix composites. J Mater Sci, 55:9848 – 9882, 2020.
[16]  S.B. Sapozhnikov and E.I. Shchurova. Voxel and finite element analysis models for ballistic impact on ceramic-polymer composite panels. Procedia Engineering, 206:182 – 187, 2017.
[17]  A. Sas, N. Ohs, E. Tanck, and G. Harry van Lenthe. Nonlinear voxel-based finite element model for strength assessment of healthy and metastatic proximal femurs. Bone Reports, 12:100263, 2020.
[18]  M.W. Beall, J. Walsh, and M.S. Shephard. Accessing CAD geometry for mesh generation. In Proceedings of the 12th International Meshing Roundtable, IMR 2003, Santa Fe, New Mexico, USA, September 14-17, 2003, pages 33 – 42, 2003.
[19]  K. Ho-Le. Finite element mesh generation methods: a review and classification. Computer-Aided Design, 20:27 – 38, 1988.
[20]  S.H. Lo. Finite element mesh generation and adaptive meshing. Progress in structural engineering and materials, 4:381 – 399, 2002.
[21]  T.R. Faisal and Y. Luo. Study of fracture risk difference in left and right femur by QCT-based FEA. Biomedical Engineering Online, 16:116, 2017.
[22]  H. Kheirollahi and Y Luo. Understanding hip fracture by QCT-based finite element modeling. J. Med. Biol. Eng, 37: 686 – 694, 2017.
[23]  Z. Yosibash, R. Padan, L. Joskowicz, and C. Milgrom. A CT-based high-order finite element analysis of the human proximal femur compared to in-vitro experiments. J Biomech Eng, 129:297 – 309, 2007.
[24]  J. Aboudi, S.M. Arnold, and B.A. Bednarcyk. Fundamentals of the Mechanics of multiphase materials. Butterworth-Heinemann Elsevier, Oxford, UK, 2013.
[25]  O. Pierard, C. Friebel, and I. Doghri. Mean-field homogenization of multi-phase thermo-elastic composites: a general framework and its validation. Composites Science and Technology, 64:1587 – 1603, 2004.
[26]  J.D. Eshelby. The determination of the elastic field of an ellipsoidal inclusion, and related problems. Proc. Royal Soc. of London, Ser. A, 241:376, 1957.
[27]  T. Gentieu, A. Catapano, J. Jumel, and J. Broughton. Computational modelling of particulate-reinforced materials up to high volume fractions: Linear elastic homogenisation. Proceedings of the Institution of Mechanical Engineers, Part L: Journal of Materials: Design and Applications, 233:1101 – 1116, 2019.