Modeling and Simulation of Fiber Orientation in Injection Molding of Polymer Composites

We review the fundamental modeling and numerical simulation for a prediction of fiber orientation during injection molding process of polymer composite. In general, the simulation of fiber orientation involves coupled analysis of flow, temperature, moving free surface, and fiber kinematics. For the governing equation of the flow, Hele-Shaw flow model along with the generalized Newtonian constitutive model has been widely used. The kinematics of a group of fibers is described in terms of the second-order fiber orientation tensor. Folgar-Tucker model and recent fiber kinematics models such as a slow orientation model are discussed. Also various closure approximations are reviewed. Therefore, the coupled numerical methods are needed due to the above complex problems. We review several well-established methods such as a finiteelement/finite-different hybrid scheme for Hele-Shaw flow model and a finite element method for a general three-dimensional flow model.


Introduction
Short fiber reinforced polymer composites are widely used in manufacturing industries due to their light weight and enhanced mechanical properties.The short fiber composite products are commonly manufactured by injection molding, compression molding, and extrusion processes.During those processes, the fibers orient themselves due to the flow and interactions between neighboring fibers and/or cavity wall.This orientation is anisotropic in general, which results in an anisotropic property of final products.Thus, the prediction of the fiber orientation during the process has been the subject of considerable amount of research during past decades.
The injection molding process is a well-established mass-production method for polymeric materials.In this process, the molten polymer or polymer composite is injected Mathematical Problems in Engineering into a mold cavity, which is namely a filling stage.After cooling the polymer material inside the mold, the final product is ejected from the mold.The overall processing time is usually less than one minute, and a complex three-dimensional shape can be produced quite easily.In the injection molding of short fiber polymer composites, the fiber orientation develops during the filling stage of the process.In this paper, we review the modeling and simulation of fiber orientation during injection molding filling stage.The prediction of fiber orientation in injection molding involves two subjects of fiber orientation simulation and injection molding simulation.In this paper, we give more attentions in the first subject of the fiber orientation while the later one is rather briefly discussed.

Fiber Orientation Kinematics
The orientation state of a group of fibers can be represented by a probability distribution function ψ p where p represents the fiber orientation vector.When the fibers are assumed to move with the bulk flow of the fluid, the conservation equation of the distribution function is written as follows 1 : where ṗ is the angular velocity vector of the fiber and D/Dt is the material time derivative.For a concentrated fiber suspension where hydrodynamic interactions and direct contacts take place between neighboring fibers, Folgar and Tucker 2 gave the following model for the angular velocity: where ω is the vorticity tensor, γ is the rate-of-strain tensor, and λ is a geometrical parameter of the particle.The last term in the right-hand side of 2.2 with D r is a rotary diffusivity term which is an additional term to the original work of Jeffery 3 to take the effect of the interactions between fibers.For an accurate calculation of the fiber orientation state, one needs to solve 2.1 and 2.2 , which requires huge computational resources for a practical implementation in injection molding simulation.For the efficiency of the computation, orientation tensors by Advani and Tucker 1 have been widely accepted, which enables a compact representation of the orientation state.The second-and the fourth-order orientation tensors are defined as follows:

2.3
The orientation tensors satisfy the full symmetric property and the normalization property as follows:

2.4
Mathematical Problems in Engineering 3 From 2.1 and 2.2 , the evolution equation of the second-order orientation tensor can be derived as follows: It should be noted that the fourth-order orientation tensor A appears in 2.5 .In a similar manner, the evolution equation of any orientation tensor contains next higher even-order tensor.Thus, one needs a closure approximation to close the set of the evolution equations of the orientation tensors.Several closure approximations are discussed in the following.

Closure Approximations
The closure approximation represents the fourth-order orientation tensor as a function of the second-order orientation.A hybrid closure approximation is a simple and stable model thus has been widely used in many numerical simulations 1 .It combines linear and quadratic closure approximations.The hybrid closure tends to overpredict the fiber orientation tensor components in comparison with the distribution function results.Also the full symmetric property of A is not satisfied due to the quadratic closure term A qua ijkl .The orthotropic closure approximations had been developed by Cintra and Tucker 4 and improved further by Chung and Kwon 5 .In the orthotropic closure, three independent components of A in the eigenspace system, namely A 1 , A 2 , and A 3 , are assumed to depend on the eigenvalues of a as follows: where a 1 and a 2 are two largest eigenvalues of a and C i K i 1, 2, . . ., 6 are eighteen fitting parameters.The orthotropic closure satisfies the full symmetry condition.There are several different versions of the orthotropic closure approximation which were developed to improve the accuracy and to overcome some nonphysical behaviors of the original model.The performance of the orthotropic closure approximation is quite successful.However, it requires additional computation for tensor transformations between the global coordinate and the principal coordinate, which is its drawback in terms of the computational efficiency.
The invariant-based closure approximation uses the general expression of the fourthorder tensor in terms of the second-order tensors of a and δ as follows 6 : where S is the symmetric operator.The six coefficients of β i are assumed to be the function of the second and third invariants of a.The invariant-based closure is as accurate as the eigenvalue-based closures, while its computational time is much less about 30% than that of the eigenvalue-based closures since it does not require any coordinate transformations.The neural-network-based closure approximation was developed recently, which assumes two-layer neural network between the fourth-and second-order orientation tensors as follows 7 : where W i are weighting coefficients, b i are biases coefficients, and f i are transfer functions.
A linear function and tangent hyperbolic function were used for the transfer function.
The neural network closure is accurate for a wide range of the flow fields, and also its computational time is much lower than the orthotropic closures.However, it requires huge numbers of coefficients for W i and b i , which can be quite troublesome for a practical application.The quantitative comparison of the computational cost and the accuracy between different closure approximations could be found in 6-8 .There were attempts to solve the evolution equation for the fourth-order orientation where one needs the closure approximation for the sixth-order orientation tensor 9 .The invariant-based approach was employed, and more accurate prediction of the fiber orientation could be achieved.However, there still remains a trade-off between the accuracy of solution and the additional computational cost to solve more equations.

Slow Orientation Model
The kinematic equation of 2.5 has been widely accepted in the numerical simulations in the past decades.However, there are some experimental observations that the actual fiber orientation kinematics might be two-to ten-times slower than the model predicts 10, 11 .One simple way to describe the slow orientation kinematics was to scale the velocity gradient with a constant factor 10, 11 .This idea is based on an assumption that the effective velocity gradient experienced by the group of fibers is less than the bulk velocity gradient due to the cluster structure of the fibers.A modified model is written as follows: where 1/κ is, namely, the strain reduction factor.However, this model is not objective, which means that it is coordinate dependent, thus cannot be used for a general flow field.
Recently, an objective model for the slow orientation was developed 12 .In this model, the growth rates of the eigenvalues of a are multiplied by κ and the resulting evolution equation is as follows: where L and M are defined in terms of eigenvalues and eigenvectors of a.The only difference with the standard model of 2.5 appears in the fourth-order orientation tensor term and the fiber interaction term.This model could fit the experimental data of transient shear viscosity, especially the shear strain where the viscosity has a peak, better than the standard model.Also the fiber orientation prediction results using the slow orientation model were in a better agreement with the experimental result 13 .Figure 1 shows the effect of slow orientation on the fiber orientation result in a center-gated disk geometry 14 .

Fiber Interaction Model
In most of the previous works, the rotary diffusivity D r of the following form has been widely used 2 : where C I is the fiber interaction coefficient and γ is the generalized shear rate.There are several models regarding C I depending on the fiber concentration.Bay 15 had carried out extensive experimental works and suggested a fitting curve for C I as a function of φ f L/D as follows: This result shows the screening effect of the fiber interaction in the concentrated suspensions.Ranganathan and Advani 16 proposed a theoretical model using Doi-Edwards theory as follows: where K is a proportionality constant and a c is the average interfiber spacing.In this model, the fiber interaction depends on the orientation states via a c .Particularly, for fiber suspension in a viscoelastic media, Ramazani et al. 17 modified 2.13 as follows: where c is the polymer conformation tensor and n is a constant.According to this model, the fiber interaction decreases as the polymers are stretched in the direction of the fiber orientation.Park and Kwon 18 also used 2.14 in developing a rheological model for a fiber suspension in a viscoelastic media, and the coupling effect between fiber and polymer in C I was found to be dominant only at the high shear rate regime.Also there had been several anisotropic diffusivity models developed in the literature 19-21 .The anisotropic diffusivity model could fit the experimental data better than the isotropic model particularly for a long fiber composite 21 .

Governing Equation
The fiber orientation develops during the filling stage of the injection molding process where the flow can be assumed to be incompressible.In addition, the inertia is negligible because of the high viscosity of the polymer melt.As a result, the mass and momentum conservation equations can be written as follows: where u is velocity vector, p is the pressure, and τ is the stress tensor.The energy conservation equation is as follows: where ρ is the density, c p is the heat capacity, η is the viscosity, and k is the heat conductivity.The last term in the right-hand side comes from the viscous dissipation.For a thin cavity geometry which is a usual situation in injection molding process, one can simplify the mass and momentum conservation equations using Hele-Shaw approximation 22 .Then, one obtains the so-called pressure equation as follows: where S is the fluidity and b is the half-gap thickness of the cavity.It might be mentioned that a generalized Newtonian constitutive relation has been used to derive 3.3 which will be discussed in Section 4. Also the energy conservation equation is rewritten as The in-plane velocities of u and v can be written in terms of the pressure gradient as The Hele-Shaw flow model have been successfully employed in the injection molding simulation and the fiber orientation prediction in the past decades 23-28 .
However, it should be noted that the three-dimensional flow at the melt front region, namely, the fountain-flow, could not be described by Hele-Shaw approximation, which can result in an inaccurate prediction of fiber orientation at the melt front region and skin layer.Figure 2 shows the effect of fountain flow in fiber orientation prediction 8 .Only a limited number of studies could be found which includes the fountain-flow effect since one should solve the three-dimensional governing equation to describe the fountain flow, which requires huge computational costs.Instead of solving the three-dimensional equations, Bay and Tucker 29 and Han and Im 30 introduced a special treatment particularly at the melt front region to describe the fountain flow.Ko and Youn 31 studied the details of gap-wise distributions of the fiber orientation for simple geometries without Hele-Shaw approximation.Han 32 and Chung and Kwon 33 studied the fountain flow effect using finite element method in an axisymmetric geometry of center-gated disk.VerWeyst et al. 34 used a hybrid approach where three-dimensional flow simulation is carried out locally with the help of Hele-Shaw flow simulation providing the boundary conditions.

Constitutive Model
A generalized Newtonian model for polymer melts has been widely accepted for injection molding simulation, which can be written as follows: This model is simple and accurate for injection molding process where the shear deformation dominates the flow 22 .There are several models for shear thinning viscosity of the polymer Mathematical Problems in Engineering melt such as the power law model, the Cross-model, and so on.The temperature dependency of the shear viscosity can be described by means of an Arrhenius model or WLF model.The generalized Newtonian model has been well adopted for fiber orientation prediction during injection molding in many literatures.Rigorously speaking, however, the fiber orientation is coupled with the momentum balance equation via constitutive relation.This coupling effect can be important depending on the flow regime 27, 28, 33, 35 .The stress tensor where the fiber anisotropic contribution is taken into account can be written as follows: where N p is a particle number.It should be noted that this model is for slender fibers where fiber thickness can be ignored and η is the viscosity of the neat polymer matrix.
There are several models for N p depending on the fiber concentration 34-37 .In practical applications in injection molding simulation, the Dinh and Armstrong model 35 could be readily employed because of simplicity in spite of its inaccuracy.The effects of coupling on the fiber orientation and flow have been studied in 26, 27, 31, 33, 38 .In the injection molding simulation, it was observed that the coupling effect is important particularly at the core and transition layers near the gate region 27, 33 .There are some models in which both the polymer viscoelasticity and fiber anisotropy are taken into account in the stress tensor 17, 18, 39, 40 .In this case, the stress tensor is written as follows: where η s is the solvent viscosity which is Newtonian , η m is the viscosity of matrix which is non-Newtonian , η p is the polymer viscosity, θ is the polymer relaxation time, and c is the polymer conformation tensor.Generally, the evolution equation of the polymer conformation tensor has the following form: where the left-hand side is the upper-convective time derivative of c and the right-hand side represents the dissipation process of the polymer.As noted in 13 , this viscoelastic model would not be useful practically in predicting the fiber orientation since already simple models of the generalized Newtonian model 4.1 or fiber suspension model 4.2 can give accurate predictions of the flow and fiber orientation results.Instead, this model would be useful to study the viscoelastic deformation or residual stress of the polymer.

Numerical Methods
The prediction of fiber orientation during injection molding involves a coupled analysis of flow, temperature, moving free surface, and fiber orientation evolution.This coupling between the equations can be treated in an explicit manner at each step of time integration or melt front advancement.Commonly, the flow field is solved first, then the temperature and fiber orientation are solved using the flow field solution.Finally, the melt front is advanced and then the next time step begins.Particularly in the momentum and energy conservation equations, there is nonlinearity due to the viscosity which depends on the shear rate γ and the temperature T .The nonlinearity has been solved by the iteration method at each time step.For the moving free surface problem, there are several different methods in the literature.When one uses Hele-Shaw flow model, the melt front advancement is simply carried out based on the nodal value of the pressure.For each node, a scalar quantity, namely, fill factor, is defined, which is updated for each time step according to the pressure solution.
For three-dimensional flow model, the free surface capturing is carried out by solving additional problem of convection equation.In this method, an artificial scalar quantity, namely, pseudo-concentration, is defined for each node so that the fluid region and the air region can be distinguished depending on the value of pseudo-concentration.

Hele-Shaw Flow
For the flow and temperature fields simulation, a finite-element/finite-difference coupled scheme has been successfully employed for injection molding simulations when the Hele-Shaw flow is assumed with a generalized Newtonian constitutive model 22 .In this method, the in-plane discretization is achieved by the finite element method while the thickness discretization is achieved by the finite difference method.For instance, the pressure and temperature fields can be approximated as follows: where N i is the finite element interpolation function of ith node, p i is the pressure at the ith node, and T l i is the temperature at the ith node of lth finite difference layer.Commonly, a triangular linear interpolation is used for N i .Since the thickness directional discretization is independent of the in-plane discretization in this method, one can have a fine discretization for the thickness directional distribution of physical quantities such as temperature, viscosity, and fiber orientation tensor components, which is quite important in the injection molding simulation.The finite difference method appears where the physical quantity depends on the in the thickness direction z of the cavity.For instance, the conduction term in the energy conservation equation when the conductivity is assumed to be constant can be approximated by the finite-element/finite-difference scheme as follows: In most previous studies, the orientation tensor is defined only at the centroid of the element of each finite difference layer.To solve the evolution of the orientation tensor 2.5 or 2.10 , the velocity gradient should be calculated at the centroid of each element from the flow field solution.For the Hele-Shaw flow, the velocity gradient tensor is as follows: where the velocity can be calculated from 3.5 .The velocity gradients in thickness direction, that is, ∂u/∂z and ∂v/∂z, can be easily defined by the finite difference method.However, one needs a special treatment to obtain the other components.It should be noted that the velocity interpolation becomes one order lower than the pressure interpolation by 3.5 .For a triangular linear interpolation of the pressure, which is the most case in practice, the velocity is constant within an element.Thus, one needs to introduce a special method to calculate the velocity gradient components, which is usually done using the velocities of the neighboring elements.Another way to tackle this problem would be to use higher-order interpolation of the pressure, which has not been reported in the literature yet.
As for the time integration of the energy conservation equation, commonly the firstorder implicit method has been widely employed, while higher integration scheme such as fourth-order Runge-Kutta method has been used to solved the fiber orientation evolution 26-28 .For the convective terms in the energy equation and fiber orientation evolution equation, one should introduce, namely, upwinding scheme to avoid numerical oscillation or wiggling of the solution.The basic concept of the upwinding scheme is to introduce different weightings for each node depending on the sign of dot product of the velocity and the outward normal vector at each node.
The finite-element/finite-difference hybrid method has been quite successful because very fine refinement is possible in the gap-wise finite difference method independently of the finite element method.This is the most important advantage of the hybrid method in the injection molding simulation because of the high shear nature of the flow inside the cavity.

Three-Dimensional Flow
The significance of the three-dimensional flow in the fiber orientation prediction has been discussed in many studies 29-33 .For the three-dimensional flow and temperature fields model, the finite element method is most widely used in the polymer processing simulation because of its flexibility for complex geometry problems 41-45 .A mixed velocity and pressure formulation has been widely used where the incompressibility is treated by a Lagrange multiplier method.One should be careful about the interpolation functions of the pressure and the velocity due to the Babuška-Brezzi condition which states that the pressure interpolation should be one order lower than the velocity interpolation.One can also find a stabilized formulation which enables equal-order interpolation of the velocity and the pressure for a convenience of finite element formulations 46 .For the three-dimensional finite element case, the numerical formulation is rather simple than for that of the finiteelement/finite-difference hybrid scheme.However, one major issue arising with the threedimensional flow model is the huge computational cost particularly for a fine discretization in the thickness direction of the cavity, which limits the practicality of the three-dimensional model 45 .
As mentioned above for convective problems such as energy conservation equation and fiber orientation evolution equation, one should use a stabilizing method, namely upwinding scheme, to get a stable and accurate solution without wiggling.In the finite element formulation, a streamline-upwind/Petrov-Galerkin SUPG method has been adopted widely in the literature 13, 33, 42 .For a convective problem of where R is the residual and X can be either temperature or orientation tensor, the SUPG formulation can be written as follows:

R, Y
e Rζ e u • ∇Y e 0, 5.5 where Y is the weighting function and ζ e is a stabilizing parameter which can depend on the element size and characteristic velocity at the element.

Conclusions
This study reviews modeling and simulation of fiber orientation during injection molding of polymer composites.The prediction of fiber orientation requires coupled analysis between the flow, temperature, free surface, and fiber orientation.The coupling can be treated in an explicit manner for each step of time integration or melt front advancement.The nonlinearity in the equations has been solved using iteration schemes.
For the flow field during injection molding filling stage, it is commonly assumed that the inertia is negligible and the flow is incompressible.Hele-Shaw approximation is adopted to simplify the flow equation, which is accurate for most part of the flow in a thin cavity geometry.For Hele-Shaw flow model, a finite-element/finite-difference hybrid scheme is widely accepted to solve the pressure equation and the energy conservation equation.Since the discretization in the thickness direction is independent of the in-plane discretization, one can achieve a fine solution in the thickness direction which is quite important in the injection molding simulation.In the Hele-Show flow model, however, the fountain flow at the melt front could not be described, which could result in a poor prediction of the fiber orientation near the melt front region particularly in the shell and skin layers.Also Hele-Shaw model is not suitable for a general three-dimensional geometry with a varying thickness.Several finite element methods have been developed to solve three-dimensional flow field in injection molding filling simulation.Because of the computational cost, however, three-dimensional flow simulation remains as one of the major problems for the accurate prediction of the fiber orientation.
The orientation state of fibers is commonly described by the orientation tensors for the efficiency of computation.Using the orientation tensor results in a need for a closure approximation to close the set of evolution equations of orientation tensors.Various closures from different bases have been developed to achieve both the accuracy and the computational efficiency.Invariant-based closure is as accurate as orthotropic closures while its computational cost is much less than orthotropic closures.Neural network closure might be as good as invariant-based closure, but it requires a lot of fitting coefficients.As for the kinematics of the fiber orientation, a modified Jeffery model by Folgar and Tucker has been a basis to understand the fiber orientation behavior.There are recent progresses in the fiber kinematic equation such as slow orientation model and anisotropic rotary diffusivity model.The slow orientation model enables more accurate prediction of the orientation state especially near the gate region.The anisotropic rotary diffusivity model would be useful particularly for the long fiber composites.The initial condition of the orientation state also can affect the final results of the simulation depending on the cavity geometry.Thus, the simulation including the sprue region could improve the accuracy of the fiber orientation prediction.
The fiber orientation can be solved either in a coupled manner or in a decoupled manner with the flow field depending on the type of constitutive model.The coupling effect is significant particularly near the core and transition layers.Most of the commercial simulation software use the decoupled analysis because of the convenience in developing each problem solver independently like a module.The coupled analysis would be more accurate than the decoupled analysis, but the constitutive model and the numerical scheme become more complicated.

11 Figure 1 :
Figure 1: Effect of slow orientation κ on fiber orientation prediction in a center-gated disk at radial position of r/b 12.4 where b is the half-gap thickness of the disk.

11 Figure 2 :
Figure 2: Effect of fountain flow on fiber orientation prediction in a center-gated disk at radial position of r/b 22.8 where b is the half-gap thickness of the disk.