Abstract
A solid rocket motor is one of the critical components of solid missiles, and its life and reliability mostly depend on the mechanical behavior of a composite solid propellant (CSP). Effective mechanical properties are critical material constants to analyze the structural integrity of propellant grain. They are estimated by a numerical method that combines the Voronoi cell finite element method (VCFEM) and the homogenization method in the present paper. The correctness of this combined method has been validated by comparing with a standard finite element method and conventional theoretical models. The effective modulus and the effective Poisson’s ratio of a CSP varying with volume fraction and component material properties are estimated. The result indicates that the variations of the volume fraction of inclusions and the properties of the matrix have obvious influences on the effective mechanical properties of a CSP. The microscopic numerical analysis method proposed in this paper can also be used to provide references for the design and the analysis of other large volume fraction composite materials.
1. Introduction
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) [1].
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. [7] 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. [8]. Another group of models published devoted to estimate mechanical properties of propellants in recent years. Yang and Liu [9] use coarse triangle meshes and ANSYS to predict the elastic modulus of the composite solid propellant. Zhi et al. [10] 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 [11] and Matous et al. [12] 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 [13] 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 [14]. In the following years, the developments of VCFEM were presented for linear elastic problems [15] and elastic-plasticity problems [16], as well as failure analysis [17]. 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 [20] and the effective properties [21] 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. [1] 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 [7]. 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 () 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 [22].(2)The number of inclusion is generated randomly from the range of to . and can be obtained using the equation: where and are the maximum value and the minimum value of mean values of inclusion radius, respectively, and 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 μm to 310 μm. Since the volume fraction is about 83%, the and 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 μ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 [14], 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: where the variables with the superscripts and correspond to the interior of the matrix and inclusion phases, respectively, while those with superscripts and refer to the variables on the element boundary and internal matrix-inclusion interface, respectively. is the stress field within the matrix or the inclusion. is the elastic compliance matrix; represents the area of the matrix or the inclusion. , , and 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. represents the prescribed tractions on the boundary . is the outward normal unit vector of the element boundary or matrix-inclusion interface.
The displacement field is interpolated by using the nodal displacements and the boundary displacement interpolation functions , 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 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 , which contains polynomial terms in the x and y coordinates variables as reported in [14], and unknown vectors of coefficients . The can be simplified as
The matrices and are defined as where
Considering the stationary condition of the total complementary energy, we can get . Consequently, the vector is expressed as
Substituting (7) into (3), with respect to the displacement , the expression of the stiffness matrix of the element is obtained as
Assembling stiffness matrices of each element,
The nodal displacements are the solutions to the following:
The Lagrange multiplier method is used to impose additional constraints to avoid rigid body displacement at the interface lMT. The inner node displacements can be represented by the element boundary nodal displacements , 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: where where is the matrix to restrain the rigid body displacements and can be expressed as
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 [7]: where () is the quantity of triangular (quadrangular) elements, () is the mean stress/strain of each single triangular (quadrangular) element, () is the area of the mth (nth) triangular (quadrilateral) integration cell, and 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 , 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 (refer to (8)). And the rank of matrix is equal to the number of the columns of matrix (refer to (5)). For elements with more nodes, it is necessary to increase to obtain enough rank. The value of is adopted 25 to analyze octagon element. When is equal to 25, the highest order of the power of matrix is equal to 4. As shown in (5), the Gauss integral with eight points is used in the calculation of matrix to conquer the ill-conditioned problem of the stiffness matrix [23].
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 Pa and Poisson’s ratio . The inclusion material is the following: Young’s modulus Pa and Poisson’s ratio . A uniform tensile stress 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 μm,0.6 μm), (1.8 μm,0.6 μm), (1.2 μm,1.2 μm), (0.6 μm,1.8 μm), and (1.8 μm,1.8 μ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).
(a) Mesh of VCFEM
(b) Mesh of FEM
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 [24]. 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 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: GPa and . The material of the inclusion is SiC: GPa and .
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.
(a) Young’s modulus
(b) Shear modulus
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 MPa and Poisson’s ratio . The AP particle material: Young’s modulus, GPa, Poisson’s ratio . The Al material: Young’s modulus GPa, Poisson’s ratio . 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 μ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 μ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 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 VAP/VAl varies with the quantity of AP cell elements xAP. The 3D graphs in Figure 7 display the relationship between the effective properties and the volume fraction of inclusions and VAP/VAl.
Figure 8(a) illustrates an upward trend with the increase of volume fraction of inclusion . Upon further analysis, we can note that the effective modulus is closed to the modulus of the matrix when the is less than 60% and, with the increase of , the effective modulus increases more rapidly. Figure 8(b) illustrates a downward trend with the increase of . The ratio of VAP/VAl describes little influence on the effective properties.
(a) Effective Young’s modulus versus Vf and VAP/VAl
(b) Effective Poisson’s ratio versus Vf and VAP/VAl
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 EM and matrix Poisson’s ratio νM.
It can be seen in Figure 9(a) that the effective modulus has a positive relation with EM, while the influence of νM can be ignored. Figure 9(b) shows that the effective Poisson’s ratio has a positive linear relation with EM. Furthermore, the effective Poisson’s ratio decreases slightly with the increase of EM. But with the increase of the matrix modulus, the decreasing trend is more and more inconspicuous.
(a) Effective Young’s modulus versus EM and νM
(b) Effective Poisson’s ratio versus EM and νM
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 EI and the inclusion Poisson’s ratio νI.
(a) Effective Young’s modulus versus EI and νI
(b) Effective Poisson’s ratio versus EI and νI
From Figure 10(a), the effective modulus increases sharply with the increase of EI when the ratio of EI/EM is less than 10. As EI continues to increase, the effective modulus keeps constant steadily. The variation of ν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 EI. The effective Poisson’s ratio increases linearly with νI. But the slope of the effective Poisson’s ratio and νI is declining with the increase of EI.
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.
4. Conclusion
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 microstructural 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.