Numerical simulation of the injection molding process of short fiber composites by an integrated particle approach

In this study, we propose an integrated particle approach based on the coupling of smoothed particle hydrodynamics (SPH) and discrete element method (DEM) to predict the injection molding process of discrete short fibers. The fibers in the coupled SPH-DEM model are treated as non-rigid bodies to allow deformation and fracture. The interaction between resin and fibers is solved by a physical model to take into consideration of drag forces. Two cases of injection molding process with different volume fractions of short fibers are studied to predict the flow behaviors of fibers and resin. The numerical results qualitatively agree with previous experimental studies. It is found that the velocity contour of resin flow is parabolic in shape due to the velocity gradient near the wall boundaries and consequently the moving direction of fibers is in parallel with the flow direction of resin. Fiber accumulation is found in the case with higher content of short fibers.


Introduction
Injection-molded short fiber composites are commonly utilized in automobile and lightweight structures, due to their distinctive advantages, such as superbly flexible molding, short time and low cost of processing, and good stiffness and strength to a certain extent [1,2]. The microstructure of such composites, e.g., orientation, length, and distribution of fibers, has significant effects on the mechanical properties of the composites. While the microstructure largely depends on the molding process, affected by the injection speed, mold geometry, and temperature. Among the effects of the microstructure, the orientation of fibers was firstly studied by introducing an orientation tensor [3], which was then widely used to predict the microstructure of the injected molding products [4][5][6]. In addition, the fiber motion can be analyzed by developing some numerical approaches, which are based on the orientation tensor and fluid dynamics. However, the fiber orientation-based approaches cannot present the microstructure of the composite in detail as the tensors cannot predict the motion of an individual fiber. Sun et al. [7] and Oumer et al. [8] simulated injection molding process using a computational fluid dynamics (CFD) code MoldFlow, and they found that in the vicinity of the mold cavity the fiber density is very large in the flow direction but very small in the thickness direction, which were validated by the experiments. Modhaffar et al. [9] simulated and characterized the fiber suspension during the material flow within the rectangular cavities with their developed CFD model, which has the ability of describing the temperature profile and predicting the fiber orientation during the injection process. Thi et al. [10] simulated the injection molding process to predict the fiber orientation in short-glass fiber composites with different fiber weight concentration. The approach was based on the Folgar and Tucker equation [11] and the Jeffrey's equation. The former one is well-known for modeling; the fiber orientation and the latter one is aiming at the interaction between fibers. These numerical results were validated quantitatively by experiments using high resolution 3D X-ray computed tomography.
To understand the mechanism of the accumulation of injected fibers within the flow, it is necessary to investigate Ke Wu and Lei Wan contributed equally to this work. the motion behavior of the fibers during the injection process. Yamamoto and Matsuoka [12,13] proposed a particle method named particle simulation method (PSM) to simulate the motion of fibers, in which all the fibers were modeled as an assembly of particles, and the equations of motion for each particle were solved separately. In their modeling, the accumulation and deformation of fibers was represented by particles, however, the resin flow and fiber motion were analyzed separately. Thus, interaction between fiber motion and resin flow as well as the influences of fiber orientation and motion on the flow field could not be considered. Yashiro et al. [1,2] utilized the moving particle semi-implicit (MPS) method [14,15] to simulate the injection process of short fiber composites. This method is capable of tracking the free surface flow, such as flow-fronts and the motion of each fiber, which is represented by an assembly of particles. However, the fibers are assumed as rigid bodies, which differs from the experimental observations where fibers could deform and even break up during the filling process [16]. In addition, both the interaction between the fibers and interaction between fibers and resin are mathematically incorporated in the MPS algorithm without a physical interpretation. For example, the interaction between fibers due to collision is handled by indirect contact due to the nature of MPS. What is more, as both resin and fibers are represented by two types of particles in MPS scheme without allowing any overlap, thus the volume fraction of fibers is constrained and there would be insufficient amount of particles to accurately capture the flow of resin when the fiber volume fraction is high.
Recently, a weakly compressible SPH-like particle method was utilized to model the flow of injection with the fiber distribution obtained and subsequently the elastic moduli of the polymerized material was computed according to the fiber orientation and concentration [17]. The approach is capable of identifying the domains where the fibers are stuck or broken and the domains where the different failures of matrix can be observed. Also, experiments were conducted for the validation of the numerical simulations. He et al. [18,19] conducted 2D and 3D SPH simulations of the injection molding process of short fiber composites, considering the polymer as a power law fluid and the fibers as rigid bodies. In the 2D simulations, a U-shape fountain flow was found in the front of the flow and the center of U-shape gradually sank during the injection process; the fibers were found to move aligning with the flow direction and accumulate at some points in the cavity. Also, the final fiber orientation was slightly influenced by the initial fiber configuration. In the 3D simulations with three layers of fibers in the thickness direction (i.e., one core layer and two skin layers), the fibers were aligned with the flow direction in the skin layer and were inclined to the flow direction in the core layer, due to the larger shear rate close to the mold wall. Though the injection molding process of short fiber composites was modeled with SPH method, the fibers were still considered as the rigid bodies, resulting in the discrepancy of the results.
Considering the research gap mentioned above, we propose a coupled SPH-DEM approach to simulate the injection process of short fiber composites. The fibers are made up of bonded DEM particles, and the resin is represented by discrete SPH particles. This coupled SPH-DEM as a meshless method is under Lagrangian scheme, thus mesh generation and remeshing can be avoided to save computational cost. The approach has been comprehensively examined and validated in our previous works [20,21] for modeling general interactions between fluids, particles, and structures. Physical models will be applied to compute the interaction between fibers and the interaction between fibers and resin. In addition, to save the computing cost the local averaging technique is implemented in this study that both the continuity equation and the interaction force, i.e., drag force, are related to the local mean voidage. Meanwhile, as deformation, collision, and/or breakage of the fibers are considered in the injection process simulations, more accurate results of fiber orientation and distribution can be obtained.
This paper is organized as follows. Section 2 explains the coupled SPH-DEM model. Section 3 evaluates two cases of the injection molding process with different content of short fibers and two cases with different configurations of mold by applying the coupled SPH-DEM model. The predicted fiber motion is qualitatively compared with the other experimental and numerical results, and the mechanism of fiber orientation and distribution is discussed. Finally, some conclusions are drawn in Section 4.

Coupled SPH-DEM model
The injection molding process of short fiber composites is a type of multiphase flow problems where short fibers are fully immersed in the resin. Modeling the entire process is very complicated, and the interaction between short fibers and resin must be considered. For simplicity and straightforward implementation in MPS modeling, each short fiber is treated as a rigid body, and the interaction between fiber and resin is calculated using the same algorithm as for the resin [1,2]. This assumption lacks of physical explanation of interaction force between fibers and resin, and to better simulate the injection molding process a coupled SPH-DEM model is proposed herein, where each short fiber is represented by four identical bonded DEM particles to allow deformation or even fracture of fibers when their ultimate bond strength is exceeded. In addition, the resin-fiber interaction is considered as fluidsolid interaction using a drag force model together with a local average technique [22].
To study the injection molding process of short fiber composites, two Lagrangian particle methods, which represent resin and short fibers respectively, are integrated. SPH treats the resin in fluid phase as a set of discrete particles moving governed by the Navier-Stokes equations. It should be noted that only particles j located in the support domain have influence on particle i as shown in Fig. 1. Four DEM particles are bonded together to represent one short fiber using a linear parallel bond model as shown in Fig. 2a, b. The springdashpot model with a frictional element is adopted to calculate the collision between fibers in both normal and tangential directions as shown in Fig. 2c. Given that each fiber is a homogeneous and well-connected granular assembly with mainly axial deformation, it can be treated by an isotropic material, which mechanical properties is described by the elastic constants of Young's modulus (E) and Poisson's ratio (ν) [23]. E and ν are emergent properties that can be related to the effective modulus (E * ) and the normal-to-shear bond stiffness ratio (k * ) at the contact as follows: E is related to E * with E increasing as E * increases, and ν is related to k * with ν increasing up to a limiting positive value as k * increases. These relationships are obtained by specifying E * and k * as follow: where R 1 and R 2 are the radius of particle pair and k * is set up to 0.25 in this study. It can be assumed that the effective modulus E * is equal to Young's modules E when the particles' stiffness is much smaller, thus the forces are predominantly carried by the parallel bonds, i.e., k n ¼ 0:01k n .
To couple the fibers and resin, the governing equations and the interaction forces are reformulated using a local volume fraction and a weighted average algorithm. More details can be referred to our previous papers [20,21], including the governing equations for SPH and DEM, the bond model with fracture criteria, the local averaging technique, and a series of validation tests. In this section, more attention will be drawn on how this coupled model is implemented to simulate the injection molding process.
The density of each SPH particle can be approximated by using a continuity density method based on the continuity equations and some transformations [24]. Note that the density can be affected by the surrounding particles within the support domain: where subscript i stands for particle i, while subscripts j stand for the surrounding ones of particle i within its support domain. m j stands for the mass of particle j, v ij stands for the relative velocity between the particle pairs. W ij stands for the kernel function and its gradient determines the contribution of these relative velocities. In this study, Wendland kernel is adopted in all simulations, i.e., where q = |r|/h. |r| stands for the distance between two particles, h represents the smoothing kernel length associated with a particle and the normalization is ensured by setting up the constant C h in the form of 7/(64πh 2 ) when considering two dimensions.
In the same way as the particle approximation of density, the moment equation in SPH form can be described as follows: where subscript r stands for the resin particles, P stands for the particle pressure, Π ij stands for the viscosity term, R ij stands for the anti-clump term and F ext stands for the external forces acting on particle i. The external forces, F ext , stand for the feedback effect from the solid particle. Due to the nature of 2D simulations in this study and the majority of injection molding machines are horizontally oriented, while ignoring the effects of gravitational and buoyancy forces. Exceptions are applications such as insert molding in which vertical machines are still used to take advantage of the gravity. Consequently, the drag force is the only dominating factor hydrodynamically affecting the DEM particle motions (for fibers) and then the change on the fluid flow field (for resin) accordingly.
The injection molding is assumed to be processed at a constant high temperature (e.g., 140°C) with negligible solidification and thus heat transfer is not taken into account in this study. The motion of each particle is calculated based on the Newton's second law considering various forces as follows: where subscript f stands for the fiber, v p stands for the DEM particle velocity, F c p stands for the sum of direct contact forces from inter-particle collisions, F b p stands for the sum of forces transferred among parallel bonds and F d p stands for the fluidsolid interaction force from the SPH side.
Direct contact and bond forces in Eq. 6 are modeled following the traditional DEM manner in which the normal force is evaluated based on the overlap between the contacting particle and the shear force is calculated in an incremental fashion. There are various interaction laws to calculate the normal force and in this study we adopt a spring-dashpot system in normal and shear directions. Details of calculating the bond force and contact force are referred to [20,21,23].
In this study, the drag force model of Ergun [25] and Wen and Yu [26] is used to calculate the fluid-solid interaction force as follows [22]: where β f is the interphase momentum transfer coefficient, v r is the average resin flow velocity around DEM particle f for fiber and then v r −v f À Á is the relative velocity between resin flow and fiber.
The value of β f can be divided into two regimes according to ϵ f : where μ f stands for the fluid viscosity, ρ f stands for the reference fluid density, C d stands for the drag coefficient of an isolated DEM particle, and d f stands for the particle diameter. The velocity of surrounding resin flow is approximated by using Shepard filter: where v r is the velocity of SPH particle. The drag coefficient C d is calculated as follows: The aforementioned drag force acts on the DEM particles and the equal and opposite forces return to the SPH particles. Employing the kernel function, these reaction forces returned to each SPH particle is determined by a partition of the drag force in proportion to the weight of each SPH particle: Replacing the external force term, F ext , in the momentum equation for fluid phase, Eq. 5 becomes: The interaction forces between particle elements in the SPH-DEM model are schematically shown in Fig. 3.
As detailed in our previous study, boundaries for SPH and DEM are treated separately. A sketch map of these treatments is given in Fig. 4. It is worthwhile stressing that DEM particle elements have no interaction with the SPH boundary though no slip boundaries are imposed on the two phases, respectively, and there may be an overlap between them.

2D simulation of injection molding process of short fiber composites
Simulations of the injection molding process of short fiber composites was conducted in two dimensions with the proposed SPH-DEM model, in which the gravitation is not considered. As mentioned above, the injection process is the critical step to determine the quality and mechanical properties of the injected products. The same injection molding process simulated by MPS in [2] is adopted for the purpose of comparisons. The configuration of the mold is depicted in Fig. 5, where the resin is considered as an incompressible viscous fluid and represented by discrete SPH particles, and the short fibers are made of four bonded DEM particles. Initially, the short fibers are generated with a random distribution within a 30mm × 30mm square resin bath in the direction vertical to the injection direction with a volume fraction of 3.8%; hence, the total quantity of short fibers is 175, same as number of the fibers used in [2]. A rigid wall at the top of the resin bath moves downwards with an initial velocity of 1 m/s to push resin and short fiber composites into a nozzle area. The velocity of the moving rigid wall is gradually decreased to ensure a stable injection molding process. Through a stable injection, the resin and short fibers pass a narrow gate in width of 6 mm and length of 5 mm to completely fill up a slim plate mold with 80 mm in width and 4 mm in length. The density of resin is 900 kg/m 3 and its dynamic viscosity is 0.1Pa • s. The density of short fibers is 2540 kg/m 3 and every individual has a dimension of 1mm × 0.25mm. The critical ultimate tensile   Table 1. In Fig. 6, the transient injection molding process at different time intervals were captured by the SPH-DEM model and compared with simulation results obtained from MPS [2]. It can be found that around 0.02 s the resin started to pass the narrow gate carrying a few fibers to fill the plate, where the fiber orientation seems random as they are driven by the resin flow. In addition, when the injection process becomes stable from 0.06 s, the orientation of the fibers starts to be consistent, forwarding to the narrow gate, which is in accordance with the simulation results from MPS. In general the SPH-DEM model predicts similar flow profile as the MPS method. In the regions where the fibers are close to the boundary walls, a velocity gradient of resin flow is noticeable due to wall effect (see Fig. 7); therefore, the orientation of fibers in these regions changes significantly and eventually the fibers tend to move in parallel with boundary walls with a lower velocity. On the contrary, the orientation of fibers in the middle of the mold initially changes slightly despite of the direct contact between fibers and gradually becomes aligned with resin flow. At time 0.0175 s, two small regions in the middle of the mold are found with low velocity; this is the result of the appearances of vortex produced by the moving wall when it squeezes the resin flow in the nozzle area. After passing the narrow gate from 0.02 to 0.025 s, the resin flow and fibers start to split into two streams to fill the slim plate. Similarly, the orientation of majority of fibers from 0.06 to 0.12 s is still aligned with the  Figure 8 shows that the velocity gradient of resin flow is very steep due to the wall effect, and the orientation of some fibers is random because of the collision with the mold walls. Besides, a small stagnation region, where the velocity of the material flow is almost zero, is captured in front of the narrow gate and in the middle of the bottom, standing the collision of the material flow. From the comparison, it could be found that the injection process is completed at time 0.12 s, but filling process still continues until 0.165 s in the MPS simulation [2]. The difference is partly due to the velocity of moving wall, which is controlled by resin pressure at every timestep in the MPS simulation. As the threshold value of resin pressure was not given in the reference, the velocity of moving wall in the current SPH-DEM simulation is estimated to linearly decrease with an deceleration of 25 m/s 2 over the first 0.02 s and then become a constant velocity of 0.15 m/s afterwards in accordance with the displacement of moving wall in each snapshots from [2] as shown in Fig. 6. In addition, a slightly irregular parabolic profile of fountain flow in SPH particles is depicted in Fig. 9, where an unrealistic free surface and voids in the resin flow appears in both directions at the flow front after flow split. This may be due to tensile instability and free surface tension. Although the flow front is satisfactorily in parabolic shape, it occupies more spaces in the plate mold leading to a faster filling process. This issue could be mitigated by considering the surface tension as an external force in the momentum equation of SPH, which is not investigated in this study.
It can be found in Fig. 10 that the tensile stress distribution of the parallel bonds in each fiber at time 0.2 s is captured. The maximum tensile stress 936.51 MPa is found, which is far below the critical tensile strength of fiber (2000 MPa). This is due to the fact that the short fibers are mainly driven by shear flow to move with the resin; the interaction force acting on the short fibers are insignificant and not enough to break the bonds. Furthermore, as the aspect ratio of fiber is only 4, the deformation of fiber cannot be noticeably observed and the tensile stress caused by tension and bending in each fiber is almost evenly distributed. It should be noted the shear stresses in each fiber are also not noticeable, and there is little deformation perpendicular to the fiber axis.
In order to investigate the effects of the fiber volume fraction on the injection process and to compare with the experiments in [1], another simulation was conducted with a higher fiber volume fraction of 8.3%, as shown in Fig. 11. There are 375 short fibers in total during the injection process. Compared to the previous case with 3.8% volume fraction, this case with a higher fiber volume fraction share the same dynamic behavior of resin flow, and the fibers close to the skin layer of mold are almost aligned with the resin flow direction  However, more accumulation of fibers is found as shown in Fig. 11. Since the resin flow velocity is larger at a distance far from the boundary wall due to the viscosity of the resin, more frequent collision between fibers occurs in comparison to the case with volume fraction of 3.8% and this partly leads to a more random orientation of short fibers in despite of resin flow, even though the fiber orientation is mainly governed by the shear rate [3]. For these fibers close to boundary walls, their velocity is relatively lower due to the boundary effect which results in more accumulation. Consequently, composite with higher fiber content is more likely to have more disrupted orientation and irregular distributions of fibers after injection. What is more, for simulation with a high fiber volume fraction of 8.3%, the content of fibers in the right path (+y) and left path (−y) after flow split is almost the same, but slightly more fibers are observed in the right path (+y) for the simulation with a low fiber volume fraction of 3.8%. There are many factors can influence the content of fibers in either paths of mold. It can be controlled by adjusting the distribution and volume fraction of short fiber at initial stage and the velocity of moving wall (e.g., constant or variable velocity) to have an impact over the entire movement of short fibers during the injection process. In Fig. 12, only the upper part of the area marked by red square in our simulation, when the injection process is finished, is extracted for the comparison with the experiments in areas A, B, and C due to the limited short fibers and relatively dense fibers gathering at the narrow gate. It should be noted that the aspect ratio of the short fiber composite is only 4 while in the experiments, the ratio is between 50 and 100.
(a) (b) Fig. 12 Qualitative comparison between a numerical prediction and b post-processed X-ray CT image [1] of fiber orientation upon the completion of filling process Fig. 13 Configuration of 2D SPH-DEM simulation of the injection molding process with corners Considering the aspect ratio, Jeffery [27] investigated the effects of the small aspect ratio on the orientation distribution of the short fibers and found that the effects are almost the same between the ratio of 4 and 20. Therefore, the aspect ratio of 4 can be utilized to capture the dynamic mechanical response of long fiber composites. It can be seen in the experiments that the orientation of most short fibers close to the wall (areas A and C) is aligned to the resin flow direction, while in the middle of the narrow part (area B), the orientation of the short fibers seems to be random distributed and the disorder is increasingly noticeable from area E to K. Because of small number of short fibers in the 2D simulation, the obtained simulated results can only be qualitatively compared with the experimental data acquired from X-ray CT. In the zoomin view of the red square region, two short fibers marked in a red circle are mostly in disorder with maximum orientation angle compared with neighboring short fibers marked in black squares. The maximum orientation angle can be caused by either collision between particles or least velocity gradient in the middle part. It can be seen that these predictions qualitatively agree with experimental results. To more quantitatively compare with experiments, a series of 2D simulations with different initial conditions of resin bath could be carried out as in [1], and this would take in account the unpredictable influence of short fibers distributions over the fiber orientation and movement (e.g., short fibers are not evenly distributed at the start of injection process). On the other hand, 3D simulations of the injection process would more accurately capture the details of fiber flow and orientation, but at a much higher computational expense.
Despite the difference in volume fraction of short fibers, the configuration of the mold was examined next to quantify the industrial applicability of the model developed for practical injection molding processes. In practice, the configuration of the mold has many characteristics, such as free form, thickness changes, ribs, cantilever snap fits, through holes, etc., which can greatly affect the distribution of fiber in the mold. In this part, the configuration of the mold based on the previous two cases (volume fraction 3.8% and volume fraction 8.3%) was further modified by adding corners Fig. 13 and branches Fig. 14 to represent the complexity in the practical injection molding processes. In both configurations of molds, the volume fraction of fiber is 3.8% and the moving rigid wall moves in the same way as in previous cases. In Fig. 15, it shows the resin flow and its velocity distribution when the mold with corners was half filled (0.06 s) and nearly filled (0.08 s). In a similar way to the plate mold, the fibers moved aligning with resin flow direction and their orientations were found to depend on the orientation of side walls. What is more, due to the sharp turn at corners, the velocity of resin flow at corners was almost zero. In those areas, the stagnation of fibers can be observed if the volume fraction of fiber is large enough. These predictions have a good agreement with the observed resin flow [28]. Figure 16 depicts the resin flow and its velocity distribution in a mold with branches at time 0.08 and 0.12 s. As expected, the fibers moved in the resin flow direction and most of their orientations are in parallel with the orientation of side walls. Due to the existence of branches, the flow separation was observed as indicated by the red arrows. The flow separation dominates as the fibers move into branches or the end of plate mold. For fibers passing through branches, the velocity directions inclined greatly, even though they eventually moved in line with the side walls.

Conclusions
In this study, an integrated SPH-DEM model is proposed to simulate injection molding process of short fiber composites. SPH is applied to model resin flow by using particle approximations of the Navier-Stokes equations and DEM is used to model the fibers through bonded particles. Local averaging technique and physical models are employed to compute the drag force when coupling SPH and DEM together. Two cases with different volume fractions of short fibers and two cases with different configurations of molds are conducted and numerical results of the resin flow velocity and fiber orientation and accumulation are analyzed and compared to other numerical and experimental results. A few conclusions are summarized as below: (1) During the injection process, fibers near the boundary walls are in a movement aligned with resin flow due to the velocity gradient of resin flow raised by the boundary effect.
(2) The voids in the fountain front of resin flow are observed due to free surface instability. Using surface tension an external force in the momentum equation of SPH could mitigate this unexpected issue. (3) The phenomenon of fiber accumulation is noticeable in both cases with volume fractions of 3.8 and 8.3%, but more significant in the second case. (4) Predicted fiber orientation from the 2D SPH-DEM simulations qualitatively agrees with 3D experimental results. However, extension of the current 2D SPH-DEM model to 3D using GPU acceleration is essential before it becomes a more practical tool for predicting and optimizing the injection process. In addition, solidification of resin after injection may also be possibly incorporated into the model so that bond formation between fibers and resin can be investigated. Upon the listed conclusions, the realistic movement of fibers in the resin flow for molding injection can be better predicted by the proposed SPH-DEM model in comparison with the conventional methods. In addition, the bond feature in DEM allows the breakage of fibers, which could influence on the resin flow, orientation and distribution of fibers and the material properties of molded part, and should be further investigated in the future. The current model was focused on model validation; thus much work needs to be done in the future to scale up the model by GPU computing so as to enable more applications in industry. Upon the implementation of GPU acceleration, the coupled model can be extended to large-scale 3D cases, then it could be used to optimize the injection points, pressure, and even the performance of the molded parts, etc., aiming to improve the product quality at a lower cost. Ideally, the model can be extended with multiphysical features to analyze temperature distribution and phase change.