The Effect of Micromechanics Models: 2D and 3D Numerical Modeling for Predicting the Mechanical Properties of PP/Alfa Short Fiber Composites

Abstract

In the present work, we propose to confront two modeling techniques for predicting the macroscopic properties of short alfa fiber-reinforced polypropylene composites. The first modeling was a micromechanical analysis using the Mori-Tanaka, Self-consistent, Diluted, Voigt, Reuss, and Neerfeld-Hill models. The second modeling was digital, using a specific finite element technique called the Projected Fiber (PF) approach. In the framework of this study, both 2D and 3D finite element analyses based on the PF approach were used. First, we proposed an inverse approach using these analytical and finite element models to predict the Young’s modulus of alfa fiber. Then, we compared the obtained results with the experiment values available in the literature. This comparison showed that the micromechanical models underestimated the alfa fiber’s Young’s modulus, while the finite element approach, PF, allowed for good framing of the experimental values. Moreover, we investigated the effect of fiber content on the predicted elastic properties of a polypropylene (PP) matrix reinforced with randomly distributed short alfa fibers. We noticed that the Diluted model was more accurate than the Mori-Tanaka and Self-consistent methods. As for the PF approach, its estimations were close to the experimental values. For example, the Young’s modulus for the PP/alfa with a 30 wt% of fiber content was underestimated with an error of 4.3%. It is shown that the 2D PF approach can provide calculated results with sufficient prediction accuracy.

1. Introduction

For composites with randomly dispersed short fibers, the well-known micromechanical homogenization models in the literature are analytical. Finite element homogenization methods, made accessible by the increasing performance of computing power, were generally reserved for composites reinforced with well-oriented fibers [1]. Finite element models are efficient and versatile for the inspection of complex structural behavior. Thereby, many researchers have worked recently on developing finite element models for short fiber composites [1,2,3,4]. Indeed, in 2021, Narendra et al. [3] used the finite element method with a micromechanics method to study the influence of fiber geometry and volume fraction on the elastic property of a studied composite. Their results showed that all the used methods behaved in a good agreement with up to a 50% volume fraction. In 2020, Hong-Bo, and Zheng-Ming Huang [5] investigated the matrix plasticity-induced nonlinear behavior of short fiber-reinforced composites based on the extended bridging model. In 2014, Kebir and Ayad [2] developed a specific finite element approach called Projected Fiber (PF). The matrix was modeled using a triangular finite element with a three-node CST (constant strain triangle) and the fiber with a two-node truss element. Their model incorporates the random distribution of short fiber reinforcements in the composite. Then, in 2012, Cunha et al. [6] also proposed a finite element model to simulate the behavior of cracks in a composite reinforced with short steel fibers. In a similar work published in 2009, Pan and Pelegri [7] presented analytical and numerical finite element tools for the analysis and design of composite materials with randomly dispersed short fibers. In 2007, Kari et al. [1] estimated, using a finite element simulation within ANSYS code, the effective elastic properties of two composite materials: the first with short fibers randomly dispersed and the second with short fibers transversely dispersed. Their numerical predictions were in between the Hashin–Strikman bounds and close to the results of self-consistent approximation. Finally, in 2005, Doghri and Tinel [8] proposed a two-step incremental formulation based on the Mori-Tanaka homogenization to predict the elastic–plastic behavior of multiphase elasto–plastic composites reinforced with non-spherical and non-aligned short fibers under cyclic and non-proportional loadings. Their numerical simulations showed that the proposed mean-field two-step homogenization approach was able to predict the effective properties of elasto–plastic matrix composites reinforced with distributed-orientation fibers with good accuracy.

On the other hand, in engineering practice, there are many situations in which we have composite material, but we do not have information about the properties of each phase component. Moreover, the properties of composite constituents may change during the material processing [9]. Therefore, to measure the local material properties, advanced experimental techniques, for example, nanoindentation [10] or micropillar compression [11], have been developed. Another interesting way to solve this problem is by combining an optimization method with a micromechanical model to create an inverse identification approach. This issue has been raised by several researchers who proposed different inverse identification approaches either for isotropic or anisotropic materials [12,13,14]. Witold Ogierman [15] used an evolutionary algorithm with a micromechanical model to compute the elastic constants of individual material phases on the basis of known properties of composite materials. He studied two composites: one reinforced with short cylindrical fibers and the other reinforced with cubic particles. He found that his inverse identification was successful only when the properties of the composite materials with at least two different volume fractions of the reinforcement were known, otherwise the identification was ambiguous. In a similar work, Burczynski and Kus [16] analyzed composites reinforced with continuous fibers. They combined finite element analysis based homogenization with an evolutionary algorithm to conduct multiscale modeling. Therefore, although it is a widely used approach, to our knowledge, there have been no attempts to use an inverse identification method to estimate the properties of alfa fibers.

In this work, we propose to identify the effective elastic properties of a thermoplastic matrix (i.e., polypropylene (PP)) reinforced with short alfa fibers using analytical homogenization models and by finite element analysis. To this end, we used the projected fiber approach as well as the most well-known micromechanical models in the literature, namely, Voigt and Reuss, Neerfeld-Hill, Mori-Tanaka, the diluted model, and the self-consistent, which we implemented using MATLAB software. The Young’s modulus of the alfa fiber was determined using a reverse approach from these models. Then, we studied the effect of the volume fraction on the Young’s modulus for the studied composites by considering the developed analytical and numerical models

2. Materials and Methods

2.1. Material and Manufacturing Process

In this work, our choice fell on a known polymer, polypropylene. It is a thermoplastic and semi-crystalline, and it is widely used in the industry. The Young’s modulus and Poisson ratio of the used polypropylene were 1.75 GPa and 0.35, respectively. No coupling agent was used to strengthen the fiber–matrix interface. The used PP granules were extruded by the APM company. As for the fiber, we used the alfa plant from the northeastern region of Morocco, more exactly around the city of Oujda. First, the received leaves were dried in the sun for three days. Then, they were mechanically crushed and, finally, sieved. The obtained fibers had a diameter of <800 μm and a length of ~2 cm. The fibers were soaked in salted water (35 g/L) at a temperature of 60 °C for 24 h. This washing with saltwater removed dirt, dust, and some of the waxes. It thus made the fibers more open to the next treatment [17]. The second step was to wash the fibers with distilled water and then put them in a 10% NaOH solution at room temperature for 24 h. Finally, the fibers were washed several times with distilled water and dried in an oven for 12 h at 105 °C. The studied composites had different weights of treated alfa fibers (10%, 20%, 30%, and 40 wt%). We used a single screw extruder operating at a speed of 150 revolutions per minute to produce the compounds. The six heating zones of the extruder cylinder were, respectively, brought to temperatures equal to 180, 185, 190, 190, 190, and 190 °C. The wire coming out of the extruder, once cooled, was crushed to obtain compounds 3–7 mm in length. The compounds were then dried in an oven for 12 h at a temperature of 105 °C. Finally, the specimens were obtained by injection as indicated in Figure 1.

2.2. Tensile Tests

The samples of the studied materials were subjected to tensile tests following the standard ASTM D638-10 with a cross-head speed of 2 mm/min. This norm suggests dumbbell-shaped samples with a central length, width, and thickness of 105, 10, and 4 mm, respectively. A clip-on extensometer with a 50 mm gauge length was used to measure strain. For each set of samples, all the results were taken as the average value and the standard deviation of five tested samples.

2.3. Micromechanical Models

In this work, we used some of the most well-known micromechanical models in the literature. In each model, the stiffness tensor of the composite was calculated through the following equations:

$$\mathrm{Voigt}: C_c=V_f{C}_f+V_0{C}_0$$

(1)

$$\mathrm{Reuss}: C_c^{-1}=V_f{C}_f^{-1}+V_0{C}_0^{-1}$$

(2)

$$\mathrm{Neerfeld}-\mathrm{Hill}: C_c=\left(C_{c,Reuss}+C_{c, Voigt}\right)\times 1/2$$

(3)

$$\mathrm{Self}-\mathrm{consistent}: C_c=C_0+{\sum }_{ph=1}^n{f}_{ph}\left(C_{ph}-C_0\right)[I+S_{ph}{C}_c^{-1}{\left(C_{ph}-C_0\right)]}^{-1}$$

(4)

$$\mathrm{Mori}-\mathrm{Tan}\mathrm{aka}: C_c=C_0\{I+\left({\sum }_{ph=1}^n{f}_{ph}{L}_{ph} \right){{\left[I+{\sum }_{ph=1}^n{f}_{ph}\left(S_{ph}-I\right)L_{ph}\right]}^{-1}\}}^{-1}$$

(5)

$$\mathrm{Diluted} \mathrm{model}: C_c=C_0+{\sum }_{ph=1}^n{f}_{ph}\left(C_{ph}-C_0\right):{\left[I+S_{ph}: C_0^{-1}:\left(C_{ph}-C_0\right)\right]}^{-1}$$

(6)

where C is the stiffness tensor (0 for matrix, f for fiber, ph for inclusion, and c for composite), V is the volume fraction (m for matrix and f for fiber), PFh is the inclusion’s volume fraction, and S_ph is the inclusion’s Eshelby tensor. The random distribution of the fibers is represented by the n number of phases. L_ph is the pseudo-tensor of localization defined by:

$$L_{ph}=-{\left[\left(C_{ph}-C_0\right)S_{ph}+C_0\right]}^{-1}\left(C_{ph}-C_0\right)$$

(7)

These homogenization approaches were based on the definition of a material representative volume element (RVE) and considered a linear elastic behavior for the matrix and fibers. Moreover, the fiber/matrix interface was considered perfect and that there was no inter-facial decohesion. The orientation was discretized into N families. Each of these N families of reinforcement had a particular orientation and were therefore considered as N different phases.

2.4. Finite Element Modeling: The Projected Fiber (PF) Approach

2.4.1. The 2D Approach

In 2014, Kebir and Ayad [2] developed an original approach: Projected Fiber. This approach is based on a special finite element procedure associated with a random distribution of fibers to calculate the elastic properties of composites reinforced with short natural fibers. It considers the fiber’s properties (i.e., geometry after injection, orientation, fiber content, Young’s modulus, Poisson’s ratio) and the mechanical properties of the matrix (i.e., Young’s modulus and Poisson’s ratio). In a 2D small fiber composite domain, the stiffness matrix of the reinforcement, represented by a unidimensional 2-node truss finite element, is projected on that of a plane 3-node triangular finite element associated with the resin space (Figure 2). The random aspect of the small fibers is represented by the corresponding positions of the truss finite elements. A local condensation of the fibers’ DOF (degrees of freedom) was considered, decreasing dimensions of the final stiffness matrix of the composite and reducing computational times.

2.4.2. The 3D Approach

In the 3D approach, the fiber was modeled as a truss element a–b, with a constant cross-section, described within a local coordinate system x (Figure 3a). Lf represents the microscopic fiber length associated with its diameter, Df (aspect ratio Lf/Df). We introduced them as morphological data as well as the fiber’s Young’s modulus, Ef, within the element formulation. A classical linear interpolation of the displacement field, u, is proposed:

u = N₁ u_a + N₂ u_b; L_f (N₁,N₂) = (1 − x,x)

(8)

where u_a and u_b are the nodal displacements of fiber a–b, described within the local coordinate system (x). In a 3D coordinate system (i.e., X, Y, and Z), the displacements used for the simulation are U, V, and W. Associated with the fiber nodes a and b, we defined them as follows:

(U V W)^t = N₁ (U_a V_a W_a)^t + N₂ (U_b V_b W_b)^t

(9)

the PF approach was extended to three-dimensional problems using a classical tetrahedron 4-node element to model the matrix space (Figure 3b). As for the random aspect of the small fibers, it was represented by the same principal as in the 2D approach described in the work of Kebir and Ayad [2].

In the Projected Fiber approach, the fiber was modeled with a linear 2-node element (1 degree of freedom/node) merged, inside a 3-node constant strain triangle (2 DOF/node) for the 2D approach, while for the 3D projected fiber approach, it was merged inside a 4-node tetrahedron element node (3 DOF/node). Then, instead of using the Eshelby tensors or Euler angles to obtain the global composite rigidity, the PF approach considered the nodal degrees of freedom vector of the fiber as a projection on that of the resin element [2]. The random aspect of the short fibers was represented by considering 91 families. Therefore, the summation over the n heterogeneities became a summation over the n families of orientations. As the orientation was completely random, all fiber volume fractions were identical and equal to the ratio of the total volume fraction to the number of families. Both the 2D and 3D approaches use the same principle for the projection and random distribution of the fibers. The stiffness of the composite is written as:

$$C_c=C_c+{\sum }_{ph=1}^n\frac{A_f}{l_f}\left(E_f+E_m\right){\alpha }_{ph}$$

(10)

where C is the stiffness tensor (m for matrix, f for fiber, ph for families, and c for composite); A_f, and l_f are, respectively, the cross-sectional area and the length of the fiber element. For a family of fibers oriented at an angle α, the expression of the orientation tensor,${\alpha }_{ph}$, is as follows:

$${\alpha }_{ph}=\left[\begin{array}{cccc}\cos (\alpha {)}^2& \cos \left(\alpha \right)\sin \left(\alpha \right)& -\cos (\alpha {)}^2& -\cos \left(\alpha \right)\sin \left(\alpha \right)\\ \cos \left(\alpha \right)\sin \left(\alpha \right)& \sin (\alpha {)}^2& -\cos \left(\alpha \right)\sin \left(\alpha \right)& -\sin (\alpha {)}^2\\ -\cos (\alpha {)}^2& -\cos \left(\alpha \right)\sin \left(\alpha \right)& \cos (\alpha {)}^2& \cos \left(\alpha \right)\sin \left(\alpha \right)\\ \cos \left(\alpha \right)\sin \left(\alpha \right)& -\sin (\alpha {)}^2& \cos \left(\alpha \right)\sin \left(\alpha \right)& \sin (\alpha {)}^2\end{array}\right]$$

(11)

In order to compare the results obtained by the PF approach to the experimental values of the tensile test, we chose the boundary conditions of the RVE so as to simulate a tensile test. To this purpose, we applied on one of the ends of the RVE a loading along the axis (ox) and on the opposite side a simple support in order to prevent any translation along the x-axis and the y-axis. A mobile support was also applied on the lower side parallel to the x-axis (Figure 4) to prevent movements along the y-axis.

3. Results and Discussion

3.1. Estimation of the Alfa Fiber’s Young’s Modulus

It is often difficult to measure experimentally the Young’s modulus of a short plant fiber. Indeed, it must have a length over 20 mm to leave at least 10 mm free to glue it on the paper support. The glue can also sometimes coat the effective length of the fiber and distort the test results [18]. On the other hand, plant fibers have a lumen and exhibit strong geometric variations. Thus, the calculation of the fiber breaking section is a difficult task. In addition, the part from the plant where the fiber is extracted and its growth conditions are also important parameters that must be taken into account when calculating the Young’s modulus of plant fibers [19]. Finally, the composite elaboration process can alter the mechanical properties of the fiber. Indeed, after our fabrication process, the crushed alfa fibers were randomly distributed in samples, and their diameters were between 40 and 150 µm and their lengths did not reach 15 mm (Figure 5). This led to an aspect ratio of approximately 13.33, which seemed to correspond to the fiber bundles and not to the single fibers. Accordingly, the experimental measurement for the fiber’s Young’s modulus was not suitable because the larger the diameter of the studied fiber, the greater the probability of the presence of a defect and the greater the probability of obtaining a premature rupture [18].

Therefore, we have considered a numerical determination of Young’s modulus of alfa fiber. Knowing the experimental Young’s modulus of our PP/alfa composites, we used a reverse approach to the analytical homogenization models (Mori-Tanaka, self-consistent, diluted, Voigt and Reuss, Neerfeld-Hill) and 2D and 3D PF approaches to determine the Young’s modulus of alfa fiber. For this purpose, we plotted the evolution of the macroscopic Young’s modulus of the PP/alfa composite (10%, 20%, and 30% by mass) as a function of the Young’s modulus of the alfa fiber for all the homogenization models. Figure 6a shows the results for the PP/alfa composite with 30 wt%. First, we can see that the Voigt and Reuss models gave a very broad framing of the Young’s modulus of the composite. Then, for the Mori-Tanaka, diluted, and self-consistent models, we noticed that the more the Young’s modulus of the fiber increased, the more its influence on the Young’s modulus of the composite decreased, while the PF 2D and 3D approaches showed a linear evolution of the Young’s modulus of the composite as a function of that of the fiber. Using each of the curves in Figure 6a, we projected the experimental value of the Young’s modulus of the PP/alfa composite (30% by mass) on the axis representing the Young’s modulus of the fiber (Figure 6b). However, as the experimental value of PP/alfa represented only an average of over 10 test pieces, the value found by the reverse approach may be far from the real modulus of the fiber. Thus, for greater precision in the calculations, this reverse approach was also applied to the results for PP/alfa at 10, 20, and 40 wt%.

Table 1. Results of the reverse approach for all of the studied homogenization methods.

PP/Alfa	Mori-Tanaka	Self-Consistent	Diluted Model	Voigt	Neerfeld-Hill	PF 2D	PF 3D
10%	16.24	15.16	17.35	7.03	11.03	24.30	35.22
20%	18.21	15.75	21.13	7.83	12.41	30.20	40.47
30%	17.23	14.21	21.74	8.03	12.78	31.40	36.57
40%	16.44	13.15	22.80	8.27	12.81	32.63	46.59
Average (GPa)	17.03	14.57	20.75	7.79	12.26	29.63	39.71

summarizes the results found for the different homogenization methods.

In the literature, the research works that used the same extraction process as us showed that bundles of alfa fiber had a Young’s modulus between 18.2 and 28.43 GPa [17]. Thus, the values found by Voigt, Neerfeld-Hill, Mori-Tanaka, and self-consistent were not in agreement with the experimental measurements. Indeed, the latest experimental study on the properties of a bundle of alfa fibers, published in 2014 by Khaldi et al. [20], showed that the Young’s modulus of this fiber was 28.43 ± 4.07 GPa and its Poisson’s ratio was 0.34. The obtained value from the PF 2D approach was the only one falling within the error range of this study.

3.2. Influence of the Volume Fraction of the Fiber on the Young’s Modulus of the Composite

To make a comparison among the homogenization models, we plotted in Figure 7 the evolution of the Young’s modulus of the PP/alfa composites as a function of the fiber content for the analytical and numerical homogenization methods. This study was conducted with an alfa fiber Young’s modulus of 28.43 GPa. The Young’s modulus of composites reinforced with an alfa fiber content of 10–40 wt% was determined using the homogenization models. Actually, the fiber content, which can reach up to 60–70% by mass in composites reinforced with long fibers, cannot exceed 40% by mass when it comes to short fibers and polymer matrices.

From Figure 7, it can be seen that the models of Voigt, Reuss, Neerfeld-Hill, and self-consistent did not give an accurate prediction of the composite PP/alfa’s Young’s modulus. For instance, the self-consistent, Voigt, and Neerfeld-Hill overestimated the Young’s modulus of the PP/alfa 30% by 34%, 140%, 54%, respectively. While the Reuss model underestimated it by 32% in comparison to the experimental results. As for the Mori-Tanaka and diluted models, their results showed good agreement with the experimental values. This is especially true at a low fiber rate.

It was noted that for a rate of fiber higher than 10% by mass, a significant difference begins to appear between the numerical and analytical results. The values estimated by the projected fiber approach followed the experimental results more closely, while the results of the Mori-Tanaka and diluted models slightly overestimated the effective modulus of elasticity of the composite. Indeed, the error between the estimation of the 2D approach and the experimental Young’s moduli of the PP/alfa composites at 10%, 20%, 30%, and 40% by mass were 2.7%, −1.9%, −4.3%, and −7.1%, respectively. Whereas, in the same order, the error was −4%, −10%, −15%, and −18% for the 3D approach. As to the diluted model, it overestimated the Young’s modulus of the PP/alfa composites at 10%, 20%, 30%, and 40 wt% by 6.5%, 6.7%, 7.9%, and 7.6%, respectively. Thus, it can be seen that even at a high rate of fiber (40% by mass), the PF approach, especially the 2D approach, can provide calculated results with sufficient prediction accuracy.

4. Conclusions

In this study, two modeling strategies, one analytical (i.e., Mori-Tanaka, self-consistent, diluted model, Voigt, Reuss, and Neerfeld-Hill) and the other numerical (i.e., projected fiber approach) were used to predict the properties of a bio-composite reinforced with a randomly distributed short alfa fiber. First, using a reverse approach, we determined the Young’s modulus of alfa fiber according to each of these homogenization models. The confrontation between the results obtained by the reverse approach and the experimental values found in the literature showed that the analytical models underestimated the Young’s modulus of the alfa fiber, while the 2D projected fiber approach estimation was in very good agreement with the experimental data. On the other hand, the experimental Young’s modulus of the PP/alfa composites (10%, 20%, 30%, and 40% by mass) was compared with those obtained by the homogenization models. Note that analytical methods, apart from the Reuss model, overestimated the effective elasticity modulus of the studied composites. Nevertheless, it was noticed that the identification of the fiber’s and composite’s Young’s moduli by the diluted model was closer to the experimental value than for the Mori-Tanaka and self-consistent models. As for the projected fiber approach, its estimations were very close to the experimental values. The 2D approach based on an optimized random technique for fiber orientation remains more accurate than the 3D one that uses a simple generation meshing of random fibers. Developments are in progress to obtain improved generation 3D mesh. As the orientation state of the fibers had a great influence on the results of the inverse identification, we suggest the use a tomography device to identify exactly the volume fraction of each orientation. Furthermore, inverse identification based on two variables (i.e., Young’s modulus and Poisson’s ratio) will make the results more accurate and represents work in progress. This approach will be reported in a future work.