Abstract
Some fiber types have a high aspect ratio and it is very difficult to predict their composites using traditional finite element (FE) modeling. In this study, an FE model was developed to predict the anisotropy of composites reinforced by short aramid fibers. Three fiber distribution types were studied as follows: perfectly aligned, normally distributed, and randomly distributed fibers. The elastic constants were obtained, and, for different alignment angles and parameters in the fiber orientation distribution function, their numerical results were compared to those of the Mori–Tanaka model. Good agreement was obtained; thus, the employed FE model is an excellent and simple method to predict the isotropy and anisotropy of a composite with high-aspect-ratio fibers. Therefore, the FE model was employed to predict the orientation distribution of a composite fiber with a nonlinear matrix. The predicted and experimental results agree well.
1. Introduction
Short-fiber-reinforced polymer composites have been widely applied in industrial fields. In addition to the low cost and easy processing, versatility is an important characteristic that is determined by the performance of the composite’s constituents, fiber volume fraction, short-fiber aspect ratio, fiber orientation, and interface between the matrix and fiber [1]. A composite with aligned fibers can be recognized as an anisotropic material; however, a composite with random fibers is isotropic and has completely different mechanical properties. Generally, the fiber orientation varies because of different flow patterns inside the mold, processing conditions, and the rheological properties of the material [2]. Therefore, it is difficult for the fibers in a composite to be completely or randomly aligned. The mechanical properties of a real composite with short fibers are always intermediate between those of completely and randomly aligned fibers.
Many analyses and numerical methods have been used to predict the elastic properties of short-fiber-reinforced composites. Among the analysis methods, the Mori-Tanaka model has been proven to produce the most accurate predictions [3]. To assess the effect of fiber orientation on the elastic properties, classical models, such as the Eshelby equivalent inclusion, Halpin–Tsai equation, self-consistent, shear lag, and bounding models and their extensions, have been applied to various materials. Mortazavian et al. [4] used experimental methods along with an analytical approach to study the tensile strength and elastic modulus of short-fiber-reinforced polymer composites. They proved that a laminate analogy and modified Tsai–Hill criteria provide satisfactory predictions. Tian et al. [5] investigated the effect of fiber orientation on the effective elastic properties of short-fiber-reinforced metal matrix composites using two-step mean-field homogenization procedures including the D-I/Reuss, M-T/Reuss, D-I/Reuss, M-T/Voigt, and D-M/V-R models. Jiang et al. [1] discussed the effects of fiber orientation distribution on the overall elastic properties of composites using Mori–Tanaka’s method. Ogierman et al. [6] studied the influence of process-induced orientation of the reinforcement on mechanical properties using Mori–Tanaka’s micromechanical model and a numerical method. Laspalas et al. [2] investigated the mechanical behavior of short-fiber-reinforced plastic composites by applying the Tandon–Weng micromechanical model considering the local fiber orientation distribution. Chen et al. [7] determined the effective moduli of composites containing misoriented fibers based on the Eshelby–Mori–Tanaka theory. Huang [8] presented an analytical approach to evaluate the orientation effects on the elastic properties of a composite containing randomly oriented fibers. Fu et al. [9] modeled fiber orientation distributions using a laminate analogy approach.
It is convenient to use analytical approaches; however, compared to a numerical method, the microstructure field of deformation in composites cannot be directly obtained. Therefore, numerical methods which predict the anisotropy have been developed for different materials. Tian et al. [5] employed finite element (FE) methods with solid elements for short fibers to predict the effects of fiber orientation on the effective elastic properties of metal matrix composites. Ogierman et al. [6] investigated the effect of fiber orientation using a macroscale FE model with different boundary conditions. Laspalas et al. [2] studied the elasticity and failure of short-fiber-reinforced composites using FE simulation including an orientation average. Jansson et al. [10] employed the FE method to predict the anisotropic and nonlinear behavior of plastics reinforced by glass fibers. However, in FE modeling with multiple inclusions, “jamming” of fibers [11] easily occurs, in which the generated new fibers intersect the existing fibers. Extensive efforts have been devoted to developing methods for generating the representative volume elements (RVEs) of composites reinforced by inclusions of different shapes. Yi et al. [12] proposed modified random sequential adsorption (RSA) algorithms to generate an RVE of a random chopped fiber-reinforced composite material with straight and curved fibers. Harper et al. [11] presented a method to generate RVEs with random discontinuous carbon fibers in the plane and no limitation to the volume fraction of the fiber. Li et al. [13, 14] presented modified RSA algorithms for generating RVEs with complex microstructures such as fibers and spatially randomly distributed pores based on the microstructure information of carbon/carbon (C/C) composites.
However, most numerical models have only focused on composites with a low-aspect-ratio fiber. The aspect ratio refers to the ratio of the fiber length to the diameter. For a high aspect ratio, particularly when the fiber radius is at a micron scale and the fiber length is at a millimeter scale, such as an aramid fiber, it is very difficult to establish a numerical FE model using a traditional method particularly when the volume fraction of the fibers is high. Only a few papers have reported on this topic [15, 16]. In our previous study [17], a numerical model was proposed to investigate the mechanical properties of an aramid-fiber-reinforced rubber composite with a high aspect ratio. It was proved that the FE model could simulate large deformation mechanical behavior well. However, the applicability of this proposed FE model in describing the anisotropy of a composite remains unknown.
In this study, the FE model was further developed to predict the anisotropy of short-fiber-reinforced composites. This paper is structured as follows: In Section 2, the analytical methods for predicting the elastic modulus of the composites with perfectly aligned, completely random, and probability orientation distribution fibers are introduced. In Section 3, the FE model and the algorithms of the different fiber orientation distributions are presented. In Section 4, the boundary conditions for the different mechanical behaviors are introduced and a means of obtaining the effective elastic constants is presented. In Section 5, the numerical results of the elastic constants with a different orientation distribution are compared to those of the Mori–Tanaka model and a high coincidence is shown. Finally, based on the aforementioned study, the FE model was applied to composites with a nonlinear matrix to predict the fiber orientation distribution. The numerical results were compared to the experimental results and found to agree well. This study verifies that the employed FE model is excellent to predict the elastic properties of a short-fiber-reinforced composite.
2. Analytical Models
2.1. Mori–Tanaka Model
In the Eshelby theory of equivalent inclusion, by replacing the elastic moduli of the inclusion with those of the matrix, an inhomogeneity problem can be equivalent to a homogeneity problem. Meanwhile, an eigenstrain is introduced into the inclusion domain. Based on the Eshelby equivalent inclusion theory [18], the Mori–Tanaka model [19] considers the mutual effect between the inclusions in the composite, thereby improving the accuracy of the prediction results. It has been proven that the Mori–Tanaka model is a reliable and simple method for inclusion problems and has been widely used to predict effective elastic moduli.
It is supposed that the fiber axes coincide with the principal direction of the composite matrix. The composite is subjected to a uniform stress on the far field boundary, which is regarded as the average composite stress. For a homogeneous material, the constitutive relation is , where is the elastic constant tensor of the composite. The average strain of the composite is as follows [19]:where and are the average strain of the matrix and inclusions, respectively. is the elastic constant tensor of the matrix. C1 is the volume fraction of the inclusions. A is expressed as follows:where is the unit tensor of four orders and S is the Eshelby tensor of four orders which is related to the elastic property of the matrix and the shape of the inclusions. The analytical expressions of S for an ellipsoid inclusion problem can be found in many previous studies [20, 21]. Thus, the effective modulus of the composite can be obtained as follows:Hence, the effective elastic modulus is an explicit function that is easy to obtain.
2.2. Analytical Model with Orientation Effects
In the previous section, the deduction was based on the assumption that the fibers are oriented along the principal matrix axes. However, in practice, because of processing conditions, it is difficult to achieve a perfect alignment and there is always a fiber orientation distribution in a composite. Therefore, to determine the real effective modulus, the effect of the fiber orientation distribution should be considered. In this section, we discuss the effective moduli for three cases: a perfectly aligned fiber whose principal axes present an angle referring to the principal matrix axes, a probability fiber orientation distribution, and completely random fibers.
Generally, the orientation vector of fiber can be expressed as , where and are the Euler angles as shown in Figure 1. x1, x2, and x3 are the three directions of the global coordinates, while , , and are the three directions of the local coordinates. In this case, the general framework in the previous section remains valid. The only change was to determine the tensor of Eq. (3) in the local coordinates referring to the global structural geometry. We denoted and as the unit vectors of the global coordinates (0, x1, x2, and x3) and local coordinates (0, , , and ), respectively. A relationship between them can be constructed as follows [1]:where the transformation matrix is as follows:Therefore, for the case of perfectly aligned fibers, we obtained a relationship between the elastic tensor of the global coordinates and of the local coordinates as follows:where are the elements of the transformation matrix . and are the elements of the corresponding elastic tensor.
In the case of a misaligned fiber, the orientation average of the elastic tensor of the composite is provided by the following integration [1]:where is the probability density function of the fiber orientation. For a planar problem, Eq. (7) can be simplified combined with a discrete algorithm [5, 22] of the orientation average as follows:where is the incremental angle, is the number of the angle incrementation, , and i is the integer within (0, ). is the probability of the fiber orientation within the range of (). is the elastic constant when . Eq. (8) is a simplified algorithm that avoids the difficulty of an integral, particularly for the complex probability density function . The precision of the result of Eq. (8) depends on because an increase in will lead to high calculation precision. Previous studies show that if is appropriately selected, high-precision calculation results for elastic constants can be achieved.
In the case of randomly distributed fibers, in each angle range () is recognized to be equal. It can be calculated by . Substituting it into Eq. (8), the elastic tensor is obtained.
3. Numerical Model Development
3.1. Method
In the traditional method of finite element analysis (FEA), the solid element is always used to simulate the mechanical behavior of the fiber. However, using this method, it is very difficult to establish a composite model with a high fiber volume fraction [12]. Therefore, in [17], an FE model was introduced for a planar problem of the elastic properties of a short-fiber-reinforced rubber composite. However, thus far, we did not know the applicability of this model for the analysis of the anisotropic properties of a composite.
The FE model was established as a square whose size is [17], as shown in Figure 2. Short fibers were generated and the linear hybrid beam element (BH21H) was employed in the ABAQUS FE software. For the matrix, plane stress continuum elements (CPS8R) were used. Using the embedded technique, which has been verified to yield quite similar results with those for the traditional method [23], the mesh of the matrix and the fibers are independent of each other, as shown in Figure 3, and are more regular, leading to more precise numerical results. Another advantage is that a model with high volume fraction and high-aspect-ratio fibers is realized.
(a) Perfectly aligned fibers (=50°)
(b) Normally distributed fibers
(c) Completely random distribution
3.2. Algorithm for Fiber Generation
In this study, planar problems with three fiber distribution types, i.e., perfectly aligned, probability function distributed, and completely random fibers, were discussed. The detailed flow chart can be referred to in the literature [17]. To avoid fiber intersection, the RSA algorithm was employed using Python language in the ABAQUS software. Three random numbers, x, y, and θ, which represent the moving distance along the x direction, moving distance along the y direction, and the rotation angle of the fiber, respectively, were used. To ensure higher accuracy of the fiber volume fraction, each fiber was generated with an equal length in the model, including the fibers around the model edge. For this reason, the ranges of x and y were set as follows: , , where is the fiber length. For different fiber orientation distributions, the FE models are shown in Figure 3 and the ranges of θ are interpreted as follows:(1)Perfectly aligned fibers: θ is the orientation angle of the fibers which is invariable for all the fibers in the model.(2)Completely random fibers: θ is a random number within (, ).(3)Misaligned fibers with a probability function distribution: the incremental angle is set where is the number of the angle incrementation. First, according to a given probability function of the fiber orientation . and are the fiber amount and orientation probability, respectively, distributed within an angle range (, ). Subsequently, fibers are generated with random angles within the range (, ) and the corresponding generated fiber amount is . i is an integer between 0 and . Therefore, this algorithm is a discretized method whose precision depends on . Decreasing the enhances the calculation precision.
3.3. Probability Density Function of the Fiber Distribution
For most real composites reinforced with aramid fibers of a fixed aspect ratio, there are always different fiber orientation distributions that can be described by a symmetric function, such as a normal function, and a nonsymmetric function, e.g., a Gamma, Weibull, or exponential function. To study the effect of fiber orientation distribution on the elastic properties of the composite, a normal distribution function was considered in this study. For the planar problem, only one variable, the Euler angle , is included such that the function can be expressed as follows [1]:where () is the mean and ( > 0) is the variance of the distribution. Figure 4 shows the fiber distribution when in Eq. (9) ; ; and 30, 100, and 500 (degree measure was used for ). The diagram shows that as increases, the fiber distribution is increasingly disordered. There are two extremes: when , the fibers perfectly align, and when is sufficiently large, the fiber distribution is completely random. Indeed, when increases to a large value, the range of , (, ) is unable to accommodate all the fibers because of the normal distribution characteristics. However, fortunately, the effect of the orientation angle of the fiber on the mechanical properties of the composite is periodic, which indicates that all the fibers of a orientation angle exert the same effect. According to this, can be flexibly calculated for large .
4. Numerical Implementation
4.1. Model Parameters
In this study, aramid fiber, which is widely applied in industrial products, was employed for the investigation. Typically, aramid fiber has a high aspect ratio because its diameter is very small at approximately 0.012 mm. For chopped fibers, the length is changed to meet product requirements. Because an aramid fiber is short and rigid, it is straight at the beginning. Therefore, a straight line can be employed to simulate its initial elastic property in the FE model. According to the product manual, a fiber of modulus = 70 GPa was selected for this study and the Poisson ratio was 0.3. For the polymer composite, a matrix modulus of 1 GPa was employed. The mechanical properties of the composites containing the aforementioned components were similar to those of a composite with a resin matrix, which is widely used in the industry. The volume fraction of a fiber was defined by . t is the model thickness, set to be 0.1, and is the total fiber amount in the model. According to [17], with increasing , the numerical results will converge. =30 was used in this study.
4.2. Boundary Conditions
The FE model in this study can be regarded as at a mesoscopic scale. If the model with a minimal size has the same effective mechanical properties as the bulk material, it can be recognized as the RVE. Generally, for a microscopic model, a periodic boundary condition will be employed [24]. In the literature [11], Harper et al. adopted periodic boundary conditions and provided effective material properties within an inner RVE region. In [17], considering the Hill condition [25], four boundary conditions were imposed, i.e., uniform traction, uniform displacement, mixed boundary condition, and periodic boundary condition, for the model uniaxial tension. It has been proven that, with increasing model size, results with different boundary condition tend to be the same.
Therefore, we believed the uniform displacement boundary condition would provide the most reasonable and correct numerical results for the elastic property of uniaxial tension. The boundary conditions for the different elastic behaviors are described in Table 1. The variable u is the displacement load of the node on the corresponding edge; the superscript represents the load direction. The subscript i or j represents any node on the corresponding edge and thus the expression implies that the displacement of each node is equal. For simple shear deformation, implies that the displacement along the x direction of the node on the edge is equal to that of the corresponding node on the edge. Among all the variables listed in Table 1, only and are provided, which are fixed values for each elastic behavior. For different elastic behaviors, the deformations are shown in Figure 5. Their elastic constants can be calculated by extracting the results using Python as follows:
(a) Uniaxial tensile deformation along the x direction
(b) Uniaxial tensile deformation along the y direction
(c) Shear deformation
(1) Elastic Moduli and where is the stress in the x direction of the model imposed by uniaxial tension. and are the strain in the x and y directions, respectively. RF is the reaction force of the model edge which is loaded on and n is the node amount of one edge, which is the same for each edge of the square model. and are the reaction force and displacement of the node on the edge, respectively. The elastic modulus and Poisson ratio in the x direction can be expressed, respectively, as follows: and .
(2) Elastic Moduli and . Replacing x with y in Eqs. (10), (11), and (12), and are obtained.
(3) Shear Modulus where and are the model shear stress and strain. Other variables have been previously defined.
5. Verification and Discussion
5.1. Perfectly Aligned Fibers
To investigate the applicability of the proposed model for predicting the anisotropic properties of composites with perfectly aligned fibers, models with a different fiber orientation angle relative to the x axis were established. The boundary conditions shown in Table 1 were imposed for the different elastic deformations and the results were compared to those of the Mori–Tanaka model. Figure 6 shows this comparison when =1 mm and =0.01. It can be seen that all the elastic constants in the numerical results, , , , , and , agree well with the results of the Mori–Tanaka model.
(a) , , and
(b) and
For the elastic modulus , assuming that the uniaxial tension along the x direction is imposed, the minimum occurs at approximately but not at . To understand this, the matrix and fiber stresses were investigated. By averaging the stress of the matrix and fiber of the RVE, the contributions from the resin and fiber were both considered as follows [11]:where , , and are the average stress of the model, beam, and matrix, respectively, and V is the RVE volume. The subscripts a and b are the number of beams and matrix elements, respectively, and the subscript i is their integration point number. SF is the section force of the beam which is multiplied by ; its component in the loading direction x can then be obtained. is the stress component of each integrated point of the matrix material in the loading direction. is the cross-sectional area of the beam element. is the integrated section volume of each beam element. IVOL is the parameter of the integration point volume of the matrix element.
For the composite with =1 mm and =0.01, , , and were obtained as shown in Figure 7. It can be seen that when , as increases, the fiber stresses decrease, which weakens the fiber reinforcement effect. Until increases to approximately , the stress in the fibers is compressive. Meanwhile, for the average matrix stress, it decreases when and then increases from nearly . Thus, the total stress minimum occurs at approximately , which explains the minimum of .
5.2. Different Fiber Orientation Distributions
To express the anisotropy of the composites, a parameter was introduced that refers to the deviation degree from isotropy. The isotropic factor can be defined as follows [26]:where is the corresponding element in the stiffness matrix and is expressed by . For the planar problem, only the factor was used to describe the deviation degree from isotropy. The closer is to 1.0, the more isotropic the composite is.
Models with different ( > 0) were established when in Eq. (9) was . implies that the fibers are perfectly aligned. If increases to a sufficiently large value, the fibers are considered to be randomly distributed. For one , in the same manner, the elastic constants, , , , , , are calculated via simulation and an analytical model. Then, is obtained. The relationship between and the five elastic constants is plotted as shown in Figure 8. The comparisons between the numerical and Mori–Tanaka model results are shown to agree well. As shown in Figure 8, the minimum longitudinal modulus occurs when is near 1 indicating a randomly distributed fiber. Meanwhile, the transverse modulus, , and shear modulus, , both show a maximum value. In addition, with decreasing , seldom changes at first but the effect of on is greater than that on .
(a) , , and
(b) and
In addition, the relationship between the isotropic factor and the is shown in Figure 9. It can be seen that when increases to approximately 2000, is approximately 1.0 and tends to remain constant, which indicates a randomly distributed fiber. Figure 9 can serve as a reference to determine the degree of the anisotropy of composites and is applicable to all short-fiber-reinforced composites of normal orientation distribution. In addition, other anisotropic models with different fiber microstructural parameters were calculated and there was good agreement as shown in Figures 6 and 7. Therefore, we believe that the employed FE model predicts the anisotropy of a composite very well.
6. Application to Nonlinear Material
6.1. Effect of Orientation Distribution
The hyperelastic property of rubber composites reinforced by aramid fibers has been studied [17]. By averaging the experimental stress–strain curves along the x-direction (along the flowing direction) and y-direction (vertical to the flowing direction) specimens, the averaged curve can be recognized as the mechanical response of a composite with randomly distributed fibers. By comparing to the experimental data, the numerical model of isotropy with random fibers predicts the nonlinear elastic response of stress–strain well. For real specimens, because the short fibers in the composite are always oriented along the compound flowing direction, anisotropy always exists; however, only by means of macromechanics, it is difficult to determine the orientation degree. As described in this section, the anisotropic FE model was extended to nonlinear material.
The material parameters of the two components, rubber and short aramid fibers, are the same as [17]. The Ogden N3 hyperelastic constitutive model was employed to describe the rubber matrix elastic behavior. Based on the tensile experiment results, the material parameters were obtained as follows: =-8.116, =1.975, =5.725, =2.388, =5.467, and =-4.495. For the aramid fibers, the elastic modulus is 132 GPa, the fiber diameter is 0.012 mm, and the studied fiber length is 1.5 mm. The difference in this study was the composite fiber distribution, which was established with an anisotropic morphology. Still, a perfect interface and straight morphology for the aramid fibers at the initial deformation were supposed. For fiber length =1.5 mm and fiber volume fraction =0.0049, using the anisotropic FE model, the stress–strain response of the composites with different perfectly aligned angles was obtained as shown in Figure 10. The curve shapes with blue scatters shown in Figures 10(a) and 10(b) are similar to those of and of the linear material shown in Figure 6, respectively. The minimum longitudinal modulus shown in Figure 11 also occurs when . The stress–strain responses of the composites with different from the numerical results are shown in Figure 11, which corresponds to Figure 8 for the linear material.
(a) Stress–strain response, when the composites were loaded along the x direction
(b) Stress–strain response, when the composites were loaded along the y direction
(a) Stress–strain response when the composites were loaded along the x direction
(b) Stress-strain response when the composites were loaded along the y direction
6.2. Determination of the Orientation Distribution
For the hyperelastic material, the constitutive behavior is always expressed in the form of strain energy W. The hyperelastic constitutive model of an isotropic material has been comprehensively studied [27, 28]. However, for an anisotropic material, many terms remain to be explored [29–31]. To describe the orientation degree of the fibers as the isotropic factor of the linear material, a variable, was introduced, where and are the strain energy of the elastic response when the tensile load is the along x and y directions, respectively. According to the numerical results, for fiber length =1.5 mm and fiber volume fraction =0.0049, the relationships between and and are shown in Figure 12. It can be seen that when is near 1.0, and coincide.
From the effect of on as shown in Figure 13, according to the experimental result , the corresponding of the fiber orientation distribution was determined. , which was approximately 1400, is the most reasonable value. Similarly, for =1.5 mm and =0.0146 was determined to be approximately 3000. The experimental results [12] and the corresponding stress–strain curves of the determined are compared in Figure 14. There are good agreements, indicating that the proposed method predicting the fiber orientation distribution is applicable.
(a) =1.5 mm, =0.0049
(b) =1.5 mm, =0.0146
7. Conclusion
In this study, an FE model with an embedded technique was employed to predict the anisotropy of polymer composites reinforced with short aramid fibers. Using this FE method, one can obtain the RVE with high aspect ratio fibers in the composites. To obtain the elastic properties of the real material, the RVEs of three types of fiber orientation distributions were established using Python language: perfectly aligned, normally distributed, and randomly distributed fibers. The five elastic properties, tensile elastic moduli and , shear elastic modulus , and Poisson ratios and , were obtained by different elastic deformations with corresponding boundary conditions. For different fiber orientation distributions, the numerical results were compared to those of the Mori–Tanaka model and found to agree well.
Based on the aforementioned conclusion, the FE model was applied to predict the fiber orientation distribution. As a nonlinear material, an isotropic factor in the form of strain energy was proposed to present the deviation degree from isotropy. According to the relationships between the isotropic factor and variance, the most reasonable variances used to describe the fiber orientation distribution in the real composite were determined. Comparison between the experimental results and numerical results of the determined variance indicated good agreements. The FE model used in this study is a simple and convenient method. The numerical prediction for a composite with a high-aspect-ratio fiber has its advantages. In addition, it can be applied to simulate the large deformation of a hyperelastic material that is difficult to realize using the traditional method because of the easy occurrence of numerical nonconvergence. The proposed FE model aids in predicting the anisotropy of short-fiber-reinforced composites and has a shorter experimental period and lower cost.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
The Natural Science Foundation Projects of the Fujian Province of China (grant number 2018J01427) and the National Natural Science Foundation of China (grant number 11372074) supported this study.