Prediction of Young’s Modulus for Injection Molded Long Fiber Reinforced Thermoplastics

: In this article, the elastic properties of long-ﬁber injection-molded thermoplastics (LFTs) are investigated by micro-mechanical approaches including the Halpin-Tsai (HT) model and the Mori-Tanaka model based on Eshelby’s equivalent inclusion (EMT). In the modeling, the elastic properties are calculated by the ﬁber content, ﬁber length, and ﬁber orientation. Several closure approximations for the fourth-order ﬁber orientation tensor are evaluated by comparing the as-calculated elastic stiffness with that from the original experimental fourth-order tensor. An empirical model was developed to correct the ﬁbers’ aspect ratio in the computation for the actual as-formed LFTs with ﬁber bundles under high ﬁber content. After the correction, the analytical predictions had good agreement with the experimental stiffness values from tensile tests on the LFTs. Our analysis shows that it is essential to incorporate the effect of the presence of ﬁber bundles to accurately predict the composite properties. This work involved the use of experimental values of ﬁber orientation and serves as the basis for computing part stiffness as a function of mold ﬁlling conditions. The work also explains why the modulus tends to level off with ﬁber concentration.


Introduction
During the last decade, an enhanced demand for lightweight materials in automotive applications has resulted in the growth of the use of thermoplastic-discontinuous fiber composites [1]. The increased growth of the use of these thermoplastic matrix composite systems is due to the combination of mechanical properties and melt processability. Long-fiber (lengths > 1 mm) thermoplastic composites (LFTs) possess significant advantages over short fiber (<1 mm) composites in terms of their mechanical properties while retaining their ability to be injection molded [2]. The goal of this research is to improve the stiffness properties predictions for injection molded LFTs. During the plasticating stage of injection molding, significant fiber attrition will occur leading to a broad fiber length distribution (FLD) [3]. Fiber orientation distribution (FOD) is another highly anisotropic feature of the final injection molded parts induced by the mold filling process [4]. Mechanical properties of LFTs are highly dependent on these microstructural variables imparted by the injection molding process [5].
The computation of the elastic stiffness for the aligned and monodispersed short fiber composites was well studied by a large range of people. Tucker et al. [6] reviewed the micromechanical models for this type of composite. By comparing the standard micromechanical models with their finite element method, the authors have shown that the Halpin-Tsai equation gives reasonable estimates for stiffness, but the best predictions come from the Mori-Tanaka model based on the Eshelby's equivalent inclusion method. A similar conclusion was reached by Hine [7]. Ingber and Papathanasiou [8] found the variant of the Halpin-Tsai model is in very good agreement with their boundary element method (numerical simulation) for their entire range of fiber volume contents and aspect ratios, although they had no experimental data for comparison. These models for aligned and monodispersed fiber composites can predict the properties of a representative volume element that can subsequently be averaged to include effects of fiber length and orientation distributions of a real injection molded material. Hine et al. [7] carried out a numerical simulation using a distribution of fiber lengths generated by the Monte-Carlo technique. Garesci et al. [5] applied the fitted probability density function to get the averaged property from each single fiber length. A relatively concise way is to replace the FLD with some sort of mean fiber length [7,9,10]. The most widely used orientation-averaging scheme is to use second and/or fourth order orientation tensors developed by Tucker and Advani [4] to average the property constants. Hine et al. [11] have shown that the results determined by the constant strain orientation averaging method (assuming the units have the same strain and average their stiffness constants) were in good agreement with their numerical simulations.
There exists very little modeling work for predicting stiffness properties on injection molded LFTs [5]. The scenario of the LFTs is different from the works that studied short fiber composite. LFTs injection molding pellets prepared by the pultrusion technique have received much attention. In particular, the fibers are in the form of aligned fiber bundles coated by thermoplastic matrix, which produces pellets with much higher fiber contents than those of more conventional 'short-fiber' compounds. During the compounding process, filamentization of fiber bundles and fracture of the resultant monofilaments into elements of a lower aspect ratio lead to the dispersion of fibers into the polymer matrix [12]. For high content fiber composites, the presence of fiber bundles seems to be unavoidable which could result in a reinforcement with a much lower aspect ratio and effective stiffness than well-dispersed fibers, consequently giving a lower and even decreasing stiffness [13]. However, the existing stiffness models assume the fibers are fully and evenly dispersed in the matrix. The predicted values keep increasing with the fiber content, which is not true within commercial fiber concentration for injection molded LFTs [13,14]. In this paper, we report on the development of an empirical model to correct the fibers' aspect ratio for the actual as-formed LFTs with fiber bundles under high fiber content. After the correction, the analytical predictions show good agreement with the experimental stiffness values from tensile tests on the LFTs for the whole fiber content range investigated. Our analysis shows that it is essential to incorporate the effect of the presence of fiber bundles to accurately predict the composite properties.

Orientation
A single rigid fiber can be represented by a unit vector, p, parallel to the fiber's long axis. The average orientation state of a fiber composite can be descried by even-ordered structural tensors [4]. In this study, the second and fourth order orientation used in the modeling are defined as: where ψ(p) is the probability distribution function for orientation, and the bracket · denotes the average quantity over a volume domain. The second order orientation tensor, A, can be measured, while the fourth order orientation tensor, A 4 , can either be obtained from experiments or estimated in terms of A using various closure approximation methods [4]. In this study, several approximation closures are implemented to calculate A 4 , and the stiffness results evaluated from these approximations are compared with that obtained from the experimentally measured A 4 .

Fiber Length
Due to the compounding process, the injection molded LFTs will end up with a very broad fiber length distribution. The actual fiber length information can be described by the experimentally obtained probability of finding a fiber with length l i given by: where, N i is the measured number of fibers with length l i . A more concise approach is to replace the FLD with a single length, normally the number or weight average length defined, respectively, as:

Elastic Properties
As mentioned in the introduction, both the Halpin-Tsai (HT) and Eshelby-Mori-Tanaka (EMT) methods are extensively studied for the unidirectional, or short fiber reinforced composites. The Halpin-Tsai method is the most widely used micromechanical model because of its ease of implementation [15,16]. The corresponding equations are derived from the self-consistent ideas of Hill [17], and the final implementation is semi-empirical in nature. Several authors have concluded that a constant strain assumption works better than the constant stress assumption, that is to say, stiffness averaging surpasses the compliance averaging in the computation [11,18]. So in this article, the HT equations are used to predict the compliance matrix, and then the stiffness matrix is obtained from the inverse of the compliance matrix. Finally, the material property is calculated by averaging the stiffness constants based on FLD and FOD.
In this study, the EMT method is also used and compared with the HT method. The EMT method is a combination of Eshelby's equivalent inclusion method and Mori-Tanaka's back stress analysis, and so, this model is valid even for large volume fraction of fibers. In particular, the equivalent inclusion method of Eshelby is applied in the computation of energy-release rate in terms of the equivalent eigenstrains defined in the fiber and crack [19]. As a result, the EMT method provides the overall stiffness of the composite weakened by fiber-end cracks.
To include the FLD effect on the elastic properties, we use the EMT or HT method to calculate the stiffness of a 'reference' unidirectional fiber composite using either the experimental fiber length probability pl(l i ): or using an average fiber length L avg (L n or L w ) evaluated from the FLD: where C ijkl is the stiffness matrix having a specific fiber aspect ratio l i /d or L avg /d, and d is the single fiber diameter.
To get the stiffness of the actual injection molded LFTs, a mean tensor averaging procedure is used [20]. Specifically, the stiffness of the calculated unidirectional fiber composite including the FLD effect is averaged over the as-formed fiber orientation state as follow: where δ ij denotes the Kronecker delta and the scalar parameters B i (i =1-5) are the invariants of the stiffness tensor of the calculated unidirectional fiber composite given as: The composite specimen can be considered as a stacking sequence of thin layers which might have a different fiber orientation state. The stiffness of each layer is calculated based on the overall FLD and the characterized orientation state using the above-mentioned method. Finally, the classical lamination theory [21] is then applied to calculate the overall effective engineering stiffness of the composite [20].

Materials and Methods
The composites under investigation were 10 wt%, 30 wt%, 40 wt%, and 50 wt% glass fibers in a polypropylene matrix. The material was received from SABIC Innovative Plastics (Ottawa, IL, USA) as 12.5 mm long pellets created through a pultrusion process in 30 wt% and 50 wt% formulations. Samples with 10 wt% fibers were diluted with neat polypropylene, while 30 wt% and 50 wt% pellets were used to create the 40 wt% composites. The pellets contain a unidirectional bundle of fibers that must be dispersed during the injection molding process, specifically in the plasticating unit. Center-gated disk geometries were formed by injection molding as shown in Figure 1. In this study, the Hele-Shaw region (60% disk radius) and the advancing front region (85% disk radius) were investigated for fiber length, fiber orientation, and mechanical properties. Tensile specimens were cut from the injection-molded samples, and the young's modulus was measured according to ASTM D3039 [22] for polymer matrix composite materials. For fiber length measurement, a method based on Reference [23] was modified as follows: Instead of directly injecting epoxy into the samples after the burning off of the polymer matrix, a needle coated with epoxy was inserted into the sample all the way through the thickness direction shown in Figure 2a. As a result, fibers at this location were collected by the epoxy on the needle as shown in Figure 2b. Due to this sampling method, there was biased toward longer fibers, and the measured fiber length was then corrected based on the length of each fiber and the diameter of the needle with epoxy coating. The needle with the fibers was re-burned to get rid of the epoxy. Loose glass fibers were then dispersed on a desktop scanner and imaged at 3200 dpi. At least 3000 fibers per sample were measured using our in-house developed MATLAB ® codes. The fiber length distribution followed the typical log-normal distribution commonly observed for fiber composites. Three samples were used to produce the averages in Table 1. The number average (L n ) and weight average (L w ) fiber lengths were calculated according to Equations (4) and (5), respectively. It can be seen that the average fiber lengths in these samples are reduced with increasing fiber content.  Measurements of fiber orientation were also made to further investigate its relationship to the stiffness performance at the same locations. Orientation measurements were taken along the r-z plane, such that r denotes the flow direction with the velocity gradient in z. Samples were polished using modified metallographic techniques and oxygen plasma etched to enhance the contrast of the glass fiber and polypropylene matrix. Details of the sample preparation and orientation measurement procedure can be found in the References [24,25]. Figure 3 shows the measured through thickness fiber orientation for various glass fiber concentrations at the Hele-Shaw region. For the fiber orientation tensor, A, the diagonal components are the most important. They describe the alignment of the population with respect to the axis of the coordinate system. A value that approaches one indicates increased alignment in that direction. Only the θ direction component is presented here, because the young's modulus was measured along this transverse direction. At 30 wt%, 40 wt%, and 50 wt%, the through thickness fiber orientation distributions are very similar showing the characteristic shell-core-shell layer structures. Generally, the θ direction component reaches its largest value near the center of the disk, because that the center of the disk is dominated by extensional flow in this tangential direction. However, at 10 wt%, the distribution deviates significantly from the rest. The through thickness fiber orientation distribution is relatively 'flat' compared to those with higher fiber concentrations. This might be due to the concentration effects on the fiber orientation dynamics. At 10 wt%, the degree of fiber-fiber interaction is much less, that is to say, hindrance to fiber alignment is much less, as a result, the tangential direction alignment is quite dominant.

Results and Discussion
In all the calculations in this article, the measured FLD or the corresponding average fiber length (L n or L w ) was input into Equations (6) or (7), while the experimental second-order tensors were used in Equation (8). To analyze the accuracy of various closure approximations for predicting the elastic properties, the fourth order orientation tensors were evaluated by the linear (LIN), quadratic (QUA), hybrid (HYB), Invariant-based optimal fitting (IBOF), and improved orthotropic (ORW3) closure approximations [4,18,26,27]. Then, the as-calculated elastic stiffness using each of these estimated fourth order tensors was compared with that using the original experimental 'true' fourth-order tensor (TRU) obtained from Equation (2). The comparisons of the effective engineering modulus along the tangential direction are presented in Figures 4 and 5 applying the methods of EMT and HT, respectively. The predicted results and the general pattern with the HT and EMT methods in the studied fiber content range are very similar. Both models show a similar linear increase in the transverse modulus with increasing fiber content. However, all the values calculated from the EMT model, no matter what length parameter and closure approximation, are slightly greater than those from the HT model. Moreover, the differences between the two methods become more notable as fiber content increases. In Tucker and Liang's [6] review of the stiffness predictions for unidirectional fiber composites, for composites with an aspect ratio larger than 10, the EMT also has predicted greater values of the dominant modulus than the HT model. To answer the question which closure method or methods are the best for stiffness prediction purposes, the results calculated using the experimental fourth order orientation tensors are used as criteria. It seems that, for the entire fiber content range and all the scenarios using different fiber length parameters (FLD, L n , and L w ), the magnitudes of IBOF and ORW3 predictions are the most comparable to the criteria. Another aspect of this paper is to examine the effects of the length parameters (FLD, L n , and L w ) on the stiffness predictions for injection molded LFTs. It is seen that, the predictions of the L n parameters lie between the largest values generated by the L w parameters and the smallest predictions from the measured FLD. This result indicates that, for the purpose of replacing the FLD by an average fiber length in the computation, the L n might outperform the L w in terms of generating a better match with the FLD.  The predictions are also compared with the transverse (θ direction) tensile test results shown in Figure 6. At most locations through the thickness, the fibers are predominantly oriented toward this direction as shown in Figure 3. Here, both HT and EMT models are applied with the measured FLD and experimental fourth order orientation tensor. Only at the lowest 10 wt% (v f = 0.0385) fiber concentration do the predictions match well with the experimental results. As concentration increases, the deviations of the predictions from the experimental data turn out to be more significant. Several authors [13,28] experimentally observed that any incremental increase in fiber content appears to bring a lower improvement in properties than the previous one. That is to say, the mechanical performance of the injection molded LFTs will reach a plateau or even decrease at very high fiber concentration range. There are several possible reasons for the degradation of the mechanical properties. First, due to the non-homogeneous nature of these materials, problems can arise during their manufacture, which result in void content/porosity in the final parts [29]. Second, there is also a possibility of poor adhesion between the glass fibers and the matrix [30]. Finally, at higher fiber content, fiber bundles are very common in the injection molded LFTs [13]. There are two major effects of the presence of fiber bundles on the mechanical performance of the composites. First, the clumping of fibers will reduce the effective fiber aspect ratios in the reinforcement. Second, fiber bundles have an effect on stress concentration. Specifically, the failed fiber will induce stress concentration in those un-failed neighbor fibers within the bundles. This stress concentration occurs during the nonlinear stage of the tensile test [19,31]. Therefore, it is valid to ignore stress concentration and exclusively consider the effects of reduced aspect ratio on modulus. In this study, the density of the 50 wt% injection molded samples was measured by the pycnometry method described in Reference [32]. The density given by the supplier is 1.33 g/mm 3 and the measured value of the injection molded center-gated-disk (CGD) is 1.327 ± 0.0175 g/mm 3 , which means the void content/porosity in the final parts is negligible. There is no information about the adhesion between the glass fibers and the matrix from the supplier. In this study, we assume the adhesion is perfect to simplify the problem, which most likely not be true. However, the 10 wt%, 30 wt%, 40 wt%, and 50 wt% materials have the same surface treatments for the fibers (they are the same series using the same formulations). Therefore, it is legitimate to only include the effects of the clumping of fibers on the level-off of the elastic properties of the injection molded LFTs as fiber content increases. The cross-sectional microscopic images on the r-z plane with glass fiber foot-prints are shown in Figure 7. It is seen that the clumping of fibers turned out to be worse as fiber concentration increased. To include the effects of reduced aspect ratio on the modulus due to the existence of fiber bundles, an empirical model was proposed to modify the effective fiber bundle diameter d c in the stiffness computation .
where, d e f f ective is the effective bundle diameter, d 0 the single fiber diameter, d c a correction coefficient, a r the fiber aspect ratio, l the fiber length, and v f the fiber volume fraction. There are two empirical parameters, which need to be determined: v c is a critical volume fraction, and n is an exponent index.
Both  This empirical model was used to correct the bundles' size and fit both HT and EMT models to the tensile test results in the Hele-Shaw region with the non-linear least squares fitting method. The measured FLD and experimental fourth order orientation tensor were also used in the calculation. The empirical parameters of the d c obtained by the fitting of both HT and EMT models are shown in Table 2. The comparisons of the fitted results with experiments at the Hele-Shaw region is shown in Figure 8a. After the application of the empirical model, the predictions turns out to be much more accurate when compared with tensile test results. However, this model might over-predict the bundles' size, because the perfect adhesion between glass fibers and matrix are assumed which might not be true. The empirical parameters obtained from the Hele-Shaw region were applied to calculate the modulus at the advancing-front. The comparisons among the as-calculated predictions, corrected predictions, and the experimental results are shown in Figure 8b. The predictions also show significant improvement after the diameter correction.  Figure 9 shows the fiber bundles' size as a function of fiber length at various fiber concentration using Equation (10). The two empirical parameters were obtained by fitting the predictions from the EMT method to experimental tensile data. At low concentration (10 wt%), the magnitude of d c barely increased in the length range from 0.06 mm to 12 mm. The slopes of the lines increase notably with concentration. At higher concentration (30,40, 50 wt%) the values of d c show a rapid increase with fiber length. The trends of d c can be explained qualitatively from the fiber breakage aspect. At higher fiber content, the dominant fiber breakage mechanisms are fiber-fiber and fiber-machine interactions [34]. The residual fiber length exhibited a linear decrease following an increase in the fiber content [14,35]. Several authors also suggested that the fiber breakage rate was proportional to the fiber length or fiber aspect ratio [12,36]. The calculated bundle size, d c , for those longer fibers is very large, especially for the 50 wt% fiber content. Qualitatively, under high fiber content, those long fibers have a great chance to contact the neighboring fibers and the machine (screw and wall of barrel). Fiber bundles reduce the effective fiber length or aspect ratio, significantly, allowing the preservation of longer fibers for higher content fiber composites during the intensive injection molding process.

Conclusions
The stiffness properties have been studied in this work for 10 wt%, 30 wt%, 40 wt%, and 50 wt% injection molded LFTs (glass/PP). Experimental measurements of fiber length distribution and fiber orientation were obtained to calculate the transverse modulus (θ direction) using both Halpin-Tsai (HT) model and the Mori-Tanaka model based on Eshelby's equivalent inclusion (EMT). It has been shown that the EMT method generates slightly larger values than does the HT model. The accuracy of the various closure approximations for predicting the elastic properties has also been evaluated. The IBOF and ORW3 approximations turned out to be the best ways of describing the fourth order tensor in terms of the second order tensor. An important finding from this work is that the models over-predicted the modulus for higher fiber content injection molded LFTs. An empirical model was developed to correct the fibers' aspect ratio in the computations obtained for the actual as-formed LFTs with fiber bundles under high fiber content. After the correction, the analytical predictions matched well with the experimental stiffness values from tensile tests on the LFTs. The leveling off of the elastic properties of the injection molded LFTs as fiber content increases (>30 wt%) is due to the existence of the fiber bundles. In this study, we assume the adhesion is perfect to simplify the problem, which most likely not be true. However, the 10 wt%, 30 wt%, 40 wt%, and 50 wt% materials have the same surface treatments for the fibers. Therefore, it is legitimate to only include the effects of the clumping of fibers for comparison purpose. Moreover, the calculated bundle sizes from the proposed empirical model can be further applied to predict the strength in future work.