Effective Mechanical Property Estimation of Composite Solid Propellants Based on VCFEM.
A composite solid propellant (CSP), a highly packed particulate composite, is a prime structural material of a solid rocket motor in addition to an energetic material. A CSP consists of polymeric binder matrix (e.g., HTPB) and particle inclusions. The particle materials are usually ammonium perchlorate (AP) and aluminum (Al) .
Parametric variations of effective material constants have been studied by numerical and analytical methods in the past [2, 3] and present [4-6]. There are a few studied on the microstructural morphology of the CSP. The effective mechanical properties of CSP are critical to study the deformation and fracture characteristics of propellant grain. Various experimental and numerical studies demonstrate that the mechanical properties of the CSP can be highly sensitive to the microstructural morphology such as the dimension, shape, distribution, and properties of the inclusion. It is a kind of natural choice to study the mechanical properties of the CSP from the aspects of microscopic mechanical. In the early stage, some classical theoretical models of conventional composites have been tailored to estimate the effective properties of the CSP. These theoretical investigations are always limited to simple geometries but are often incapable of reflecting the realistic microstructure of the CSP. Instead, the numerical method becomes increasingly popular because its analysis is based on the realistic geometric simulation model. Zhang et al.  employed a homogenization theory and the displacement finite element method (FEM) to compute the mean temperature and heat flux of the CSP's representative volume element (RVE). The effect of orientation and shape of oxidizer particles on the burning rate was examined by a direct numerical simulation approach developed by Plaud et al. . Another group of models published devoted to estimate mechanical properties of propellants in recent years. Yang and Liu  use coarse triangle meshes and ANSYS to predict the elastic modulus of the composite solid propellant. Zhi et al.  use ABAQUS to study the effects of the critical contact stress, initial contact stiffness, and contact failure distance on the damaged interface model. Matous and Geubelle  and Matous et al.  develop a damage analysis tool at multiple scales from particle packing to failure use of a numerical framework.
However, the conventional finite element method requires complex grid element and huge computational costs, which limits replication and application in microstructure analysis; furthermore, very small elements may occur owing to the fact that the space among the inclusion is too narrow to create a perfect mesh. For example, they have to shrink particles in contact with other particles or reduce the volume fraction to create high-quality meshes [7, 12]. To overcome some of the limitations discussed above, Ghosh and his coworkers proposed a new numerical method known as the Voronoi cell finite element model (VCFEM) to analyze heterogeneous materials. In 1990, they  proposed a two-dimensional automatic mesh generation technique to discrete the composite domain to yield an aggregate of convex Voronoi polygons. An assumed stress hybrid formulation has been implemented to utilize the resulting Voronoi polygons as elements in a finite element model in 1993 . In the following years, the developments of VCFEM were presented for linear elastic problems  and elastic-plasticity problems , as well as failure analysis . Over the last few years, the VCFEM was further developed to address some engineering problems [18, 19]. In addition, some contributions had been devoted by other researchers to develop this method to analyze the thermomechanical  and the effective properties  of heterogeneous materials. However, very few works have been attempted to tailor the VCFEM to estimate effective mechanical properties of the CSP. Shen et al.  introduced a noninclusion VCFEM to analyze material viscoelastic constants.
The feasibility study on the application of VCFEM in the estimation of effective mechanical properties of the CSP is carried out in this paper. We focus on establishing a numerical model for the analysis of composites with high particle volume fraction. A displacement comparison between the results of VCFEM and those of commercial FEM software is carried out to indicate the correctness of the program code. And the validity of the proposed combined method is obtained by employing two classical theory methods from the literature. In addition, a simple case is analyzed to understand the influence of microstructural morphology on the effective modulus and effective Poisson ratio of the CSP.
2. Computational Procedure
2.1. Microstructure Model of CSP. The CSP is one of the highly particulate composite materials, always with particle volume fraction between 70% and 80%, or even higher . Figure 1 is a SEM image of the HTPB propellant. Spherical particles are bonded together tightly by the polymer matrix. It is necessary to model a RVE to reflect their microstructure features before obtaining effective material constants of the CSP numerically. Voronoi tessellation is a simple but effective geometric representation for charactering the microstructures of the composite. The tessellation can subdivide a plane into many Voronoi cells determined by a set of centers. Each cell may be perceived as the intersection of open half-planes bounded by the perpendicular bisectors of lines joining one inclusion center with each of its neighbors. However, the inclusion will be dissected into multiple segments by element boundaries of the mesh, when the radius of a circular inclusion in a Voronoi cell is greater than the minimum distance from its center to the center of other cells. The centers generated randomly may lead the model to be inhomogeneous usually. If an element center is closed to one adjacent center than any other adjacent centers, the room between the two centers cannot be filled with inclusion. It will result in the failure to get a large volume fraction of circular inclusion.
In fact, the coordinates of circular inclusions are almost determined when the volume fraction of circular inclusions is large and their sizes are almost equal. A new center-generated project is proposed here. In this project, the locations of centers are generated on the basis of the corresponding existing optimal equal circle packing schemes when the center's number is determined. Those circles will be homogeneous in a RVE model. And the distances between any center and other adjacent centers are almost equal.
The procedure to build the microstructure model with a large volume fraction of inclusions is described below.
(1) The size of RVE ([L.sub.RVE]) is determined based on the inclusion size and material property. RVE sizes are defined as the minimum size of a microstructural model that meets the requirement of statistical homogeneity .
(2) The number of inclusion is generated randomly from the range of [n.sub.min] to [n.sub.max]. [n.sub.min] and [n.sub.max] can be obtained using the equation:
[n.sub.min] = [V.sub.f][L.sup.2.sub.RVE]/[pi][r.sup.2.sub.max], [n.sub.max] = [V.sub.f][L.sup.2.sub.RVE]/[pi][r.sup.2.sub.min], (7)
where [r.sub.max] and [r.sub.min] are the maximum value and the minimum value of mean values of inclusion radius, respectively, and [V.sub.f] is the volume fraction of inclusion.
(3) The radius and center positions of the circle inclusion are read from the corresponding existing optimal equal circle packing schemes. The RVE is meshed with the convex Voronoi polygons according to the centers of the system.
(4) Transform the circles into ellipses. The random number generator is used to generate the lengths of major axes, ratio of major to minor axis length, and orientation (angle) of the major axis of each inclusion. The nonoverlapping constraint has to be imposed. To establish this, each inclusion is contained inside a circle having a diameter equal to the length of the major axis.
Figure 2 describes an example. The length and width are 6 mm. The mean values of the radius of inclusions are ranging from 300 [micro]m to 310 [micro]m. Since the volume fraction is about 83%, the [n.sub.min] and [n.sub.max] are 95 and 101, respectively. We set the number of inclusion to 100. And the RVE model is meshed following the second step and third step. The radius of inclusions is 308 [micro]m, and the volume fraction is 80%. Further, the inclusions are transformed into ellipses randomly following the fourth step. We can get a model with a volume fraction of 61%.
2.2. Hybrid Element Formulation. In this section, the finite element formulation with Voronoi polygons will be reviewed briefly. Figure 3 shows an example of a hybrid element, used in the VCFEM method, with an embedded inclusion. Each node of inclusion may be perceived as the intersection of the circular inclusion and the line joining the inclusion center with each matrix node. The element formulation, as reported in , is based upon the stationary of total complementary energy principle. The total complementary energy of the element with an embedded heterogeneity is the addition of the energy of the matrix and inclusion:
[mathematical expression not reproducible], (2)
where the variables with the superscripts M and I correspond to the interior of the matrix and inclusion phases, respectively, while those with superscripts E and MI refer to the variables on the element boundary and internal matrix-inclusion interface, respectively. [sigma] is the stress field within the matrix or the inclusion. S is the elastic compliance matrix; A represents the area of the matrix or the inclusion. [l.sub.E], [l.sub.MI], and [l.sub.T] represent the displacement boundary of the element, the interface of matrix-inclusion, and traction boundary, respectively. [??] represents the compatible displacement fields on the boundary of element or inclusion. [bar.T] represents the prescribed tractions on the boundary [l.sub.T]. n is the outward normal unit vector of the element boundary or matrix-inclusion interface.
The displacement field u is interpolated by using the nodal displacements q and the boundary displacement interpolation functions L(u = Lq), while the stress components within the element are assumed to be compatible with the prescribed boundary tractions and satisfy the equilibrium conditions into the region neglecting the body forces. The stress field 0 is expressed as the polynomial functions of the x and y coordinates, by using complete forms of the stress airy functions. This results in the product of an interpolation matrix P, which contains polynomial terms in the x and y coordinates variables as reported in , and unknown vectors of coefficients [beta]([[sigma].sup.M/I] = [P.sup.M/I][[beta].sup.M/I]). The [[PI].sub.mc] can be simplified as
[[PI].sub.mc] ([beta], q) = 1/2 [[beta].sub.T]H[beta] - [[beta].sup.T]Gq. (3)
The matrices H and G are defined as
[mathematical expression not reproducible], (4)
[mathematical expression not reproducible], (5)
[mathematical expression not reproducible]. (6)
Considering the stationary condition of the total complementary energy, we can get [partial derivative][[PI].sub.mc]/[partial derivative][beta] = 0. Consequently, the vector [beta] is expressed as
[beta] = [H.sup.-1]Gq. (7)
Substituting (7) into (3), with respect to the displacement q, the expression of the stiffness matrix of the element is obtained as
[K.sub.e] = [G.sup.T][H.sup.-1]G. (8)
Assembling stiffness matrices of each element,
K = [n.summation over (i=1)] [([K.sub.e]).sub.i]. (9)
The nodal displacements are the solutions to the following:
Kq = Q. (10)
The Lagrange multiplier method is used to impose additional constraints to avoid rigid body displacement at the interface [l.sub.MT]. The inner node displacements [q.sup.I] can be represented by the element boundary nodal displacements [q.sup.E], while the inner node displacements are not affected by other elements directly. Therefore, the stiffness matrix of elements and corresponding mechanical load vectors can be eliminated to reduce the computing scale as follows:
[K.sup.*.sub.E] = [K.sub.11] - [K.sub.12][K.sup.-1.sub.22][K.sup.T.sub.22], [Q.sup.*] = [F.sub.out] - [K.sub.12][K.sup.-1.sub.22] [F.sub.in], (11)
[mathematical expression not reproducible], (12)
where [phi] is the matrix to restrain the rigid body displacements and can be expressed as
[mathematical expression not reproducible]. (13)
2.3. Homogenization Methods. As shown in Figure 3, the inclusion and matrix domains of a 2D Voronoi cell element can be divided into triangular and quadrilateral integration regions, respectively. The mean stress and the mean strain are related to the values of every integration cell and can be obtained from the homogenization theory :
[mathematical expression not reproducible], (14)
where [N.sub.tri] (Nquad) is the quantity of triangular (quadrangular) elements, [[bar.[sigma]].sup.tri.sub.m]/[[bar.[epsilon]].sup.tri.sub.m] ([[bar.[sigma]].sup.quad.sub.n]/ [[bar.[sigma]].sup.quad.sub.n]) is the mean stress/strain of each single triangular (quadrangular) element, [S.sup.tri.sub.m] ([S.sup.quad.sub.n]) is the area of the mth (nth) triangular (quadrilateral) integration cell, and [S.sub.RVE] is the total area of the RVE.
2.4. Numerical Procedures
Remark 1. Octagons are used to simulate the circular inclusions. To construct the matrix H, an integration subdivision scheme is needed to achieve the numerical area integration in (5). An integration subdivision scheme is proposed to reduce the number of integration regions and enhance the precision of numerical integration. As can be seen in Figure 3, the inclusion phase and the matrix phase is subdivided into 3 and 8 quadrilateral regions, respectively.
Remark 2. The rank of the stiffness matrix is determined by the rank of matrix H (refer to (8)). And the rank of matrix H is equal to the number of the columns of matrix P[n.sub.P] (refer to (5)). For elements with more nodes, it is necessary to increase [n.sub.P] to obtain enough rank. The value of [n.sub.P] is adopted 25 to analyze octagon element. When [n.sub.P] is equal to 25, the highest order of the power of matrix P is equal to 4. As shown in (5), the Gauss integral with eight points is used in the calculation of matrix H to conquer the ill-conditioned problem of the stiffness matrix .
3. Numerical Examples
3.1. Result Validation
3.1.1. Displacement Analysis. In order to verify the applicability of the Voronoi cell element method developed, the displacement of a heterogeneous material with five inclusions subjected to uniformly tensile load is considered in this example (Figure 4). It is computed using VCFEM, and the results are compared with those obtained from the displacement-based FEM software Nastran. The matrix material is the following: Young's modulus E = 1000 Pa and Poisson's ratio v = 0.2. The inclusion material is the following: Young's modulus E = 3000 Pa and Poisson's ratio v = 0.2. A uniform tensile stress q =10 Pa at the top of the model is considered. Displacement constraints are imposed on the bottom side of the plate in the vertical direction. The center coordinates of five inclusions are the following: (0.6 [micro]m, 0.6 [micro]m), (1.8 [micro]m, 0.6 [micro]m), (1.2 [micro]m, 1.2 [micro]m), (0.6 [micro]m, 1.8 [micro]m), and (1.8 [micro]m, 1.8 [micro]m).
The RVE model is subdivided into 5 cells used in Voronoi tessellation as shown in Figure 5(a). The RVE model is meshed by FEM in the commerce procedure Patran as well. There are 15,062 quadrilateral elements as shown in Figure 5(b).
Table 1 shows the comparisons of the displacements of nodes on matrix boundary and matrix-inclusion interface in the loading direction. The results computed by the VCFEM and Nastran show a good agreement with each other. The maximum relative error is less than 5.0%. In consideration, the stiffness matrix obtained by the displacement finite element method based on the minimum principle of potential energy is greater than the real stiffness matrix. So its displacement results are smaller than real displacement results. VCFEM is based on the minimum principle of residual energy. The stiffness matrix obtained by VCFEM is smaller than the real stiffness matrix. So its displacement results are greater than real displacement results. The error between the two displacement results is acceptable.
3.1.2. Modulus Analysis. Furthermore, another example is adopted to verify the correctness of the combined method. The effective Young's modulus and shear modulus are analyzed for several volume fractions of the aluminum-SiC composites using two theoretical models and the proposed method . The VCFEM models were generated by using the program introduced in Section 2 and consisted of 100 hybrid elements containing polygonal inclusions as shown in Figure 2. A uniform tensile stress q =10 MPa at the top of the model is considered. Appropriate displacement constraints are imposed on the bottom side of the plate in the vertical direction. The material of the matrix is aluminum: E = 70.576 GPa and v = 0.33. The material of the inclusion is SiC: E = 450.0 GPa and v = 0.17.
It is clear from Figure 6 that those curves, respectively, obtained by these three methods have the same change trend. Eshelby's model and Mori-Tanaka's model are based on simplifying assumptions. They can only determine a region of real solution. And the result obtained by VCFEM is between the results of two models. It indicates that the method can obtain an effective result and shows its potential in the prediction of mechanical properties of solid propellants.
3.2. Effective Mechanical Properties of CSP. In the following case, the effective mechanical properties of the CSP are analyzed for different volume fractions of inclusion and component material properties with the RVE model of Section 3.1.2. The matrix will be assumed as an elastic material: Young's modulus [E.sub.M] =10 MPa and Poisson's ratio [v.sub.M] = 0.495. The AP particle material: Young's modulus, [E.sub.AP] = 32.4 GPa, Poisson's ratio [v.sub.AP] = 0.14. The Al material: Young's modulus [E.sub.Al]] = 68.3 GPa, Poisson's ratio [v.sub.Al] = 0.33. The volume fraction of AP and Al particles, the modulus and Poisson's ratio of the matrix and inclusions, will be changed below to understand their relationship with effective mechanical properties. As shown in Figure 7, the effective moduli of different RVEs are calculated here to verify that the RVE size used is large enough. Keep the effective radius of particles (100 [micro]m) the same; increase the number of particles to obtain RVEs with different sizes. The results showed that when the RVE size is greater than 1500 [micro]m (contains 48 inclusion particles when the volume fraction is 60%), its effective modulus remained stable. The model used below contains 100 inclusions (as shown in Figure 2), whose size is large enough to ignore the size effect of RVE here. A uniform tensile stress q =10 MPa at the top of the model is considered. Appropriate displacement constraints are imposed on the bottom side of the plate in the vertical direction.
3.2.1. Effect of Inclusions' Volume Fraction. The mechanical properties of composites with different volume fractions of AP and Al particles are calculated to study the effect of particle volume fraction. The changing of volume fraction depends on the changing of the size of RVE. The value of [V.sub.AP]/[V.sub.Al] varies with the quantity of AP cell elements [x.sub.AP]. The 3D graphs in Figure 7 display the relationship between the effective properties and the volume fraction of inclusions [V.sub.f] and [V.sub.AP]/[V.sub.Al].
Figure 8(a) illustrates an upward trend with the increase of volume fraction of inclusion [V.sub.f]. Upon further analysis, we can note that the effective modulus is closed to the modulus of the matrix when the [V.sub.f] is less than 60% and, with the increase of [V.sub.f], the effective modulus increases more rapidly. Figure 8(b) illustrates a downward trend with the increase of [V.sub.f]. The ratio of [V.sub.AP]/[V.sub.Al] describes little influence on the effective properties.
With the increase of the volume fraction of inclusion, the mechanical properties of composites, such as a composite solid propellant, is closing to properties of the inclusion materials. When the differences between the properties of two inclusion materials are smaller than the differences between theirs and the matrix's, the variation of the proportion of different inclusions is not obvious.
3.2.2. Effect of Matrix's Material Properties. The mechanical properties of composites with different modulus and Poisson's ratio of the matrix are calculated to study the effect of the matrix material. The particle volume fraction remains at 65% invariability in this case. The 3D graphs in Figure 8 display the relationship between the effective properties and the matrix modulus [E.sub.M] and matrix Poisson's ratio [v.sub.M].
It can be seen in Figure 9(a) that the effective modulus has a positive relation with [E.sub.M], while the influence of [v.sub.M] can be ignored. Figure 9(b) shows that the effective Poisson's ratio has a positive linear relation with [E.sub.M]. Furthermore, the effective Poisson's ratio decreases slightly with the increase of [E.sub.M]. But with the increase of the matrix modulus, the decreasing trend is more and more inconspicuous.
The variation of the modulus of the matrix has an effect on both the effective modulus and the effective Poisson's ratio of the CSP, when the variation of the Poisson's ratio of the matrix has a linear effect on the effective Poisson's ratio of the CSP mainly.
3.2.3. Effect of Inclusion's Material Properties. The mechanical properties of composites with different modulus and Poisson's ratio of inclusion are calculated to study the effect of the inclusion material. The particle volume fraction remains at 65% invariability in this case. The 3D graphs in Figure 10 display the relationship between the effective properties and the inclusion modulus [E.sub.I] and the inclusion Poisson's ratio [v.sub.I].
From Figure 10(a), the effective modulus increases sharply with the increase of [E.sub.I] when the ratio of [E.sub.I]/[E.sub.M] is less than 10. As [E.sub.I] continues to increase, the effective modulus keeps constant steadily. The variation of [v.sub.I] has no significant effect on effective modulus. From Figure 10(b), we can see clearly that a dramatical change has taken place in the growth process of [E.sub.I]. The effective Poisson's ratio increases linearly with [v.sub.I]. But the slope of the effective Poisson's ratio and [v.sub.I] is declining with the increase of [E.sub.I].
Only when the modulus of inclusion is matched by the modulus of the matrix, the effective modulus and effective Poisson's ratio are affected by the modulus of inclusion significantly. However, the modulus of inclusion of the CSP is larger than the matrix's. Hence, the variation of modulus of inclusion does not have as obvious influence as expected.
The effective modulus and Poisson's ratio of the CSP, which are closely related to the volume fraction and material property of each component, are the critical material parameters to analyze the structural integrity of propellant grains. A strategy for constructing RVE models of highly packed particulate composites is presented here, which is adopted by VCFEM appropriately. A numerical programming method combined with the VCFEM and homogenization method is proposed to investigate the relationship between the micro-structural morphology and the effective properties of the CSP. Based on the examples mentioned in the above sections, the following conclusions can be drawn:
(1) The mechanical properties of the CSP are significantly affected by the volume fraction of inclusions, with the increase of the volume fraction of inclusion; the mechanical properties of composites, such as a composite solid propellant, are closing to the properties of the inclusion materials. However, the variation of the proportion of different inclusions has a minor influence. Since the properties different between the inclusion and the matrix are very large, a small change of inclusion's properties does not have a significant effect on the overall effective properties.
(2) Except that the modulus and Poisson's ratio of the matrix directly influence the effective modulus and Poisson's ratio of the CSP, respectively, the variation of the matrix modulus has modest influences on the effective Poisson's ratio of the CSP.
(3) The effective properties are affected by the modulus of inclusion significantly only when the moduli of the inclusion and matrix are close. However, as for the CSP, when the modulus of the inclusion is much larger than that of the matrix, the effects of inclusion's material properties are not obvious as expected.
Conflicts of Interest
The authors declare that they have no conflicts of interest. Acknowledgments
This study was supported by the Science Project of the National University of Defense Technology.
 L.-L. Shen, Z.-B. Shen, H.-Y. Li, and Z.-Y. Zhang, "A Voronoi cell finite element method for estimating effective mechanical properties of composite solid propellants," Journal of Mechanical Science and Technology, vol. 31, no. 11, pp. 5377-5385, 2017.
 Z. Hashin and S. Shtrikman, "A variational approach to the theory of the elastic behaviour of multiphase materials," Journal of the Mechanics and Physics of Solids, vol. 11, no. 2, pp. 127-140, 1963.
 T. Steinkopff and M. Sautter, "Simulating the elasto-plastic behavior of multiphase materials by advanced finite element techniques part II: simulation of the deformation behavior of Ag-Ni composites," Computational Materials Science, vol. 4, no. 1, pp. 15-22, 1995.
 A. Malyarenko and M. Ostoja-Starzewski, "A random field formulation of Hooke's law in all elasticity classes," Journal of Elasticity, vol. 127, no. 2, pp. 269-302, 2017.
 A. Malyarenko and M. Ostoja-Starzewski, "Spectral expansions of homogeneous and isotropic tensor-valued random fields," Zeitschrift fur angewandte Mathematik und Physik, vol. 67, no. 3, 2016.
 J. Szafran and M. Kaminski, "Bridges for pedestrians with random parameters using the stochastic finite elements analysis," International Journal of Applied Mechanics and Engineering, vol. 22, no. 1, pp. 175-197, 2017.
 J. W. Zhang, S. J. Zhi, and B. Sun, "Estimation of thermophysical properties of solid propellants based on particle packing model," Science China Technological Sciences, vol. 56, no. 12, pp. 3055-3069, 2013.
 M. Plaud, S. Gallier, and M. Morel, "Simulations of heterogeneous propellant combustion : effect of particle orientation and shape," Proceedings of the Combustion Institute, vol. 35, no. 2, pp. 2447-2454, 2015.
 J. Yang and C. Liu, "A numerical approach for the simulation of composite solid propellant," in 2009 Second International Conference on Information and Computing Science, pp. 50-52, Manchester, UK, 2009.
 S. J. Zhi, B. Sun, and J. W. Zhang, "Multiscale modeling of heterogeneous propellants from particle packing to grain failure using a surface-based cohesive approach," Acta Mechanica Sinica, vol. 28, no. 3, pp. 746-759, 2012.
 K. Matous and P. H. Geubelle, "Multiscale modelling of particle debonding in reinforced elastomers subjected to finite deformations," International Journal for Numerical Methods in Engineering, vol. 65, no. 2, pp. 190-223, 2006.
 K. Matous, H. Inglis, X. Gu, D. Rypl, T. Jackson, and P. Geubelle, "Multiscale modeling of solid propellants: from particle packing to failure," Composites Science and Technology, vol. 67, no. 7-8, pp. 1694-1708, 2007.
 S. Ghosh and S. N. Mukhopadhyay, "A two-dimensional automatic mesh generator for finite element analysis for random composites," Computers & Structures, vol. 41, no. 2, pp. 245-256, 1991.
 S. Ghosh and S. N. Mukhopadhyay, "A material based finite element analysis of heterogeneous media involving Dirichlet tessellations," Computer Methods in Applied Mechanics and Engineering, vol. 104, no. 2, pp. 211-247, 1993.
 S. Ghosh, K. Lee, and S. Moorthy, "Multiple scale analysis of heterogeneous elastic structures using homogenization theory and Voronoi cell finite element method," International Journal of Solids and Structures, vol. 32, no. 1, pp. 27-62, 1995.
 S. Ghosh and S. Moorthy, "Elastic-plastic analysis of arbitrary heterogeneous materials with the Voronoi cell finite element method," Computer Methods in Applied Mechanics and Engineering, vol. 121, no. 1-4, pp. 373-409, 1995.
 S. Ghosh and S. Moorthy, "Particle fracture simulation in non-uniform microstructures of metal-matrix composites," Acta Materialia, vol. 46, no. 3, pp. 965-982, 1998.
 S. Ghosh, J. Bai, and D. Paquet, "Homogenization-based continuum plasticity-damage model for ductile failure of materials containing heterogeneities," Journal of the Mechanics and Physics of Solids, vol. 57, no. 7, pp. 1017-1044, 2009.
 S. Ghosh, "Adaptive hierarchical-concurrent multiscale modeling of ductile failure in heterogeneous metallic materials," JOM, vol. 67, no. 1, pp. 129-142, 2015.
 L. Bruno, F. M. Furgiuele, and C. Maletta, "A hybrid method for the thermo-mechanical analysis of elastic cracks in two-dimensional heterogeneous materials," Finite Elements in Analysis and Design, vol. 43, no. 6-7, pp. 444-452, 2007.
 M. Grujicic and Y. Zhang, "Determination of effective elastic properties of functionally graded materials using Voronoi cell finite element method," Materials Science and Engineering: A, vol. 251, no. 1-2, pp. 64-76, 1998.
 M. Jiang, K. Alzebdeh, I. Jasiuk, and M. Ostoja-Starzewski, "Scale and boundary conditions effects in elastic properties of random composites," Acta Mechanica, vol. 148, no. 1-4, pp. 63-78, 2001.
 S. Moorthy and S. Ghosh, "A model for analysis of arbitrary composite and porous microstructures with Voronoi cell finite elements," International Journal for Numerical Methods in Engineering, vol. 39, no. 14, pp. 2363-2398, 1996.
 X. Liu and J. Zhang, "The study of finite element method to predict effective elastic modulus of composites," Journal of Shenyang Institute of Engineering (Natural Science), vol. 5, no. 2, pp. 175-179, 2009.
Liu-Lei Shen, Zhi-Bin Shen (iD), Yan Xie, and Hai-Yang Li (iD)
College of Aeronautics and Astronautics, National University of Defense Technology, Changsha, Hunan, China
Correspondence should be addressed to Zhi-Bin Shen; firstname.lastname@example.org
Received 20 November 2017; Revised 22 January 2018; Accepted 11 February 2018; Published 12 April 2018
Academic Editor: Filippo Berto
Caption: Figure 1: SEM image of the HTPB propellant.
Caption: Figure 2: A RVE with 100 inclusions.
Caption: Figure 3: A typical Voronoi cell element with an octagonal inclusion.
Caption: Figure 4: RVE model with five inclusions.
Caption: Figure 5: Mesh of two numerical methods.
Caption: Figure 6: Effective modulus of aluminum--SiC composites versus the [V.sub.f] of SiC.
Caption: Figure 7: Effect of RVE size on effective modulus.
Caption: Figure 8: Effective properties of solid propellants versus [V.sub.f] and [V.sub.AP]/[V.sub.AL].
Caption: Figure 9: Effective properties of solid propellants versus [E.sub.M] and [v.sub.M].
Caption: Figure 10: Effective properties of solid propellants versus [E.sub.I] and [v.sub.I].
Table 1: Comparison of two methods' displacement results in load direction. Node VCFEM ([micro]m) FEM ([micro]m) Error * Matrix node 1 0 0 -- 2 0.005820 0.005598 3.97% 3 0.011190 0.010841 3.22% 4 0.016790 0.016027 4.76% 5 0.022840 0.021831 4.62% 6 0.021620 0.020608 4.91% 7 0.021960 0.020928 4.93% 8 0.021620 0.020608 4.91% 9 0.022840 0.021831 4.62% Matrix-inclusion interface node 1 0.01584 0.01568 1.03% 2 0.01514 0.01472 2.89% 3 0.01594 0.01539 3.57% 4 0.01662 0.01634 1.72% 5 0.00543 0.00519 4.56% 6 0.00442 0.00426 3.88% 7 0.00493 0.00504 2.17% 8 0.00592 0.00591 0.14% 9 0.01058 0.01017 4.00% 10 0.00980 0.00937 4.54% 11 0.01063 0.01017 4.49% 12 0.01130 0.01100 2.73% * Error = [absolute value of VCFEM - FEM]/FEM x 100%.
|Printer friendly Cite/link Email Feedback|
|Title Annotation:||Research Article; Voronoi cell finite element method|
|Author:||Shen, Liu-Lei; Shen, Zhi-Bin; Xie, Yan; Li, Hai-Yang|
|Publication:||International Journal of Aerospace Engineering|
|Date:||Jan 1, 2018|
|Previous Article:||Long-Term Microgravity Effects on Isometric Handgrip and Precision Pinch Force with Visual and Proprioceptive Feedback.|
|Next Article:||Fatigue Delamination Crack Growth in GFRP Composite Laminates: Mathematical Modelling and FE Simulation.|