Effect of the packing stage on fiber orientation for injection molding simulation of fiber-reinforced composites

Abstract

During the packing or post-filling stage, a significant flow of polymer inside the cavity may result due to the compressibility of the polymer melt under the higher packing pressure of the injection molding process. In the meantime, the effect of the packing stage on the shell–core structure of the fiber orientation for the fiber-reinforced composites has always been a concern. Even though certain commercial packages have undergone unified simulations of the filling and packing stages, fiber orientation has usually been determined at the end of the filling stage. A recently proposed mathematic model, Improved Anisotropic Rotary Diffusion and Retarding Principal Rate, having incorporated the state-of-the-art technology of 3D injection molding simulation, has demonstrated its ability to provide reliable predictions of fiber orientation. The present numerical results concentrate on comparing and analyzing the difference in fiber orientation between the filling and packing stages, while the important effects of packing time and packing pressure are further revealed. A qualitative comparison of core thickness widths in related experimental investigations is discussed herein.

Introduction

Fiber-reinforced composites are invariably selected as familiar lightweight automotive materials due to their excellent safety and good mechanical properties. The well-known shell–core structure ¹ of fiber orientation is found in the injection molding process. Before manufacturing a part, the accurate prediction of the fiber orientation state is a critical consideration in the process of structural analysis.

The famous Folgar–Tucker equation has become a standard calculation tool to predict fiber orientation in injection-molded, fiber-reinforced parts. ¹ Several studies have shown that the predictions of the Folgar–Tucker model are in good agreement with actual experiments. In the past, many researchers used “the midplane mesh model” (also called the 2.5D or Shell model) based on the Hele-Shaw approximation employed to predict fiber orientation in injection molding simulation. The injection molding process involves both the filling and packing stages. However, the previous Hele-Shaw calculations ^{2 –4} of fiber orientation are limited to the filling stage.

In practice, the experimental data on fiber orientation are measured after the end of the packing. With reference to earlier experimental investigations, Bailey and Rzepka’s ⁵ study found that the core thickness of injection-molded long fiber composites is reduced during packing. They further explained that such an effect was due to the high viscosity gradient. So far, the effect of packing on fiber orientation has not been analyzed theoretically. Even through the use of certain commercial software, reported fiber orientation has usually been obtained at the end of the filling stage.

Generating a 2.5D mesh for a complex 3D geometric part is a very time-consuming task; in addition, the 2.5D calculation based on the Hele-Shaw approximation lacks numerical resolution. The components of a complex geometric part include various rib sizes as well as changes in the flow direction and thickness. It is quite a challenge to determine fiber orientation distributions in the requisite complex designs.

Based on the classic fiber orientation models, ^{2 –4,6 –10} Tseng et al. ^11–14 have recently developed a new objective model, Improved Anisotropic Rotary Diffusion and Retarding Principal Rate (iARD-RPR), which has been incorporated into the commercial software of 3D injection molding simulation, Moldex3D (CoreTech System Co., Taiwan). The iARD-RPR model includes only three physical parameters to enable fiber orientation predictions. Foss et al. ¹⁵ utilized Moldex3D to predict short glass fiber orientation in the end-gated plaque and the center-gated disk. Just one set of the optimal iARD-RPR parameters could yield good orientation predictions with respect to flow factors, such as fill rate and part geometry.

In the present study, we therefore use the state-of-the-art technology of 3D injection molding simulation with the objective iARD-RPR model to investigate changes in fiber orientation for a comparison of the filling and packing stages. In addition, we present an explanation of the high viscosity gradient during packing. Various packing times and packing pressures influencing fiber orientation are also analyzed herein.

Theoretical background

The flow-induced fiber orientation indicates that changes in the fiber orientation are strongly driven by the flow field. For injection-molded fiber composites, the actual flow of the fibers dispersed in the polymer melt is transient, non-Newtonian and non-isothermal, with frozen layers building up as the complex mixture flows through the mold cavity. The process is highly nonlinear, as the material properties are dependent upon the rheological and thermal conditions. For completeness, here we summarize the theoretical equations, involving fluid mechanics and fiber orientation kinetics.

Basic governing equations

The governing equations of the fluid mechanics ¹⁶ which describe the transient and non-isothermal flow motion are as follows:

$\frac{\partial \rho }{\partial t}+\nabla \cdot \rho u=0$

1

$\frac{\partial }{\partial t}(\rho u)+\nabla \cdot (\rho uu-\sigma )=\rho g$

2

$\sigma =-pI+\eta (\nabla u+\nabla u^T)$

3

$\rho C_P\left(\frac{\partial T}{\partial t}+u\cdot \nabla T\right)=\nabla \cdot (k\nabla T)+\eta {\dot{\gamma }}^2$

4

where ρ is the density, u the velocity vector, t the time, σ the total stress tensor, g the acceleration vector of gravity, P the pressure, η the viscosity, C_p the specific heat, T the temperature, k the thermal conductivity, and $\dot{\gamma }$ the shear rate.

Tait’s Pressure, Volume, and Temperature (PVT) model is to express a thermodynamic state relationship, namely, the volume of the material as a function of the temperature and pressure, ¹⁷ as follows:

$\stackrel{⌢}{V}={\stackrel{⌢}{V}}_0\left[1-C\text{ln}\left(1+P/B\right)\right]+{\stackrel{⌢}{V}}_t$

5

${\stackrel{⌢}{V}}_0=\{\begin{array}{l}{b}_{1S}+b_{2S}\overline{T},\hspace{1em}if\hspace{1em}T\le T_t\\ b_{1L}+b_{2L}\overline{T},\hspace{1em}if\hspace{1em}T>T_t\end{array}$

6

$B=\{\begin{array}{l}{b}_{3S}\text{exp}\left(-b_{4S}\overline{T}\right),\hspace{1em}if\hspace{1em}T\le T_t\\ b_{3L}\text{exp}\left(-b_{4L}\overline{T}\right),\hspace{1em}if\hspace{1em}T>T_t\end{array}$

7

${\stackrel{⌢}{V}}_t=\{\begin{array}{l}{b}_7\text{exp}\left(b_8\overline{T}-b_{9}P\right),\hspace{1em}if\hspace{1em}T\le T_t\\ 0,\mathrm{ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}}if\hspace{1em}T>T_t\end{array}$

8

$\overline{T}=T-b_5$

9

$T_t=b_5+b_{6}P$

10

$C=0.0894$

11

where these parameters, b _1S , b _1L , b _2S , b _2L , b _3S , b _3L , b _4S , b _4L , b ₅, b ₆, b ₇, b ₈, and b ₉ are fitted by experimental data.

It is important to understand the viscosity flow curve for highly concentrated fiber composites. For pure polymer fluids, the flow curve of the viscosity versus shear rate involves two regions: the Newtonian plateau at low shear rates (zero-shear-rate viscosity independent of the shear rate) and shear thinning at high shear rates (viscosity decreases with increasing shear rates); this reflects the prior understanding of general polymer rheology. Reviewing previous experimental results ^{18 –21} on fiber-filled suspensions, the researchers of rheology have demonstrated that the suspension viscosity increases with fiber content and length. For the 30 wt% long glass fiber–filled polypropylene composites with various fiber lengths of 10–30 mm, Thomasset et al. ²¹ showed that the Newtonian plateau particularly vanishes at low shear rates, instead of the yield stress viscosity, namely, a very sensitive increment of viscosity with decreasing shear rates.

In general, the Cross–William–Landel–Ferry (Cross–WLF) model ^17,22 can describe complex viscosity behaviors of polymer melts, including that the viscosity varies with shear rate for the Cross and depends on temperature and pressure for the WLF. In the earlier study of suspension rheological, the yield stress was clearly found to decrease with increasing temperature. ²³ In the present study, we therefore adopted the Yield–Cross–WLF model involves the yield stress viscosity term and the Cross–WLF model, called as the Herschel–Bulkley model ^{24, 25} :

$\eta (\dot{\gamma },T,P)=\frac{{\tau }_Y}{\dot{\gamma }}+\frac{{\eta }_0(T,P)}{1+{\left(\frac{{\eta }_0\dot{\gamma }}{{\tau }^{*}}\right)}^{1-n}}$

12

${\eta }_0(T,P)=D_1\text{exp}\left(\frac{-A_1(T-T_c)}{A_2+(T-T_c)}\right)$

13

$T_c=D_2+D_{3}P$

14

$A_2={\tilde{A}}_2+D_{3}P$

15

${\tau }_Y(T)={\tau }_{Y0}\text{exp}\left(\frac{T_Y}{T}\right)$

where seven parameters are fit by related experimental data, including n, τ*, A₁, ${\tilde{A}}_2$ , D₁, D₂, and D₃; τ_Y0 is the reference yield stress and T_Y is a characteristic temperature.

Fiber orientation equations

A single fiber is regarded as a rigid cylindrical rod. The fiber’s unit vector p along its axis direction can describe the fiber orientation. For a concise representation of the orientation of a large population of fibers, Advani and Tucker ²⁶ defined the second-order orientation tensor as:

$A={\displaystyle \oint \psi (p)pp\text{d}p}=\left[\begin{array}{ccc}{A}_{11}& A_{12}& A_{13}\\ A_{12}& A_{22}& A_{23}\\ A_{13}& A_{23}& A_{33}\end{array}\right]$

16

where ψ(p) is the probability density distribution function over the orientation space; A is the symmetric matrix and its trace is A ₁₁ + A ₂₂ + A ₃₃ = 1. Physically, A = I/3 represents the isotropic orientation state, wherein I is the identity matrix. Three diagonal components, A ₁₁, A ₂₂, and A ₃₃, correspond to the flow, cross-flow, and thickness directions, respectively.

A time evolution equation of the second-order orientation tensor is fixed on the material derivative, denoted as $\dot{A}$ . Over the last three decades, theoretical researchers in the fiber suspension rheological field have made considerable effort to determine the dynamic fiber orientation states involving short and long fibers. The modern models, based on the classic fiber orientation models, the pioneering Jeffery hydrodynamic model ⁶ and the famous Folgar–Tucker Isotropic Rotary Diffusion model, ⁷ include the Phelps–Tucker Anisotropic Rotary Diffusion model, ⁸ the Wang–Tucker Reduced Strain Closure model, ⁹ and the ARD–RSC model. ⁸ These famous and modern models are available in commercial injection molding simulation software. ²⁷

From the aforementioned fiber orientation models, Tseng et al. ^{11 –13} developed a new model, the iARD-RPR. However, the iARD tensor does not pass the classic rheological rule of Euclidean objectivity, namely, material frame indifference. It is hard to ignore the non-objective effect due to the fact that different coordinate systems can yield different answers. More recently, Tseng et al. ¹³ have improved such a non-objective issue and then proposed the objective iARD model. Here, we introduce the objective iARD-RPR model divided into three parts.

Firstly, the iARD-RPR equation contains three terms: the Jeffery Hydrodynamics ${\dot{A}}^{\text{HD}}$ , the iARD ${\dot{A}}^{\text{iARD}}$ , and the RPR ${\dot{A}}^{\text{RPR}}$ , presented as:

$\dot{A}={\dot{A}}^{\text{HD}}+{\dot{A}}^{\text{iARD}}(C_I,C_M)+{\dot{A}}^{\text{RPR}}(\alpha )$

17

${\dot{A}}^{\text{HD}}=(W\cdot A-A\cdot W)+\xi (D\cdot A+A\cdot D-2A_4:D)$

18

where ${\dot{A}}^{\text{iARD}}$ has two available parameters: the fiber–fiber interaction parameter C_I and the fiber–matrix interaction parameter C_M ; ${\dot{A}}^{\text{RPR}}$ has one parameter α, which is meant to slow down the response rate of the fiber orientation; W is the vorticity tensor; D is the rate-of-deformation tensor; $\dot{\gamma }$ is the shear rate of D, $\dot{\gamma }=\sqrt{2D:D}$ ; $\xi =(a_r^2-1)/(a_r^2+1)$ is the particle shape factor; and a_r is the fiber aspect ratio. Since a_r >> 1, ξ is assumed to signify unity. The fourth-order orientation tensor A₄ is determined by using a higher order polynomial closure approximation in terms of the second-order orientation tensor A, such as the eigenvalue-based optimal fitting closure ²⁸ and the invariant-based optimal fitting closure. ²⁹ Note that the accuracy of IBOF is as good as EBOF, and IBOF requires less computational time to obtain a solution. ²⁹

Secondly, it is significant that the rotary diffusion tensor D _r depends on the square of the objective rate-of-deformation tensor for defining a new iARD model, ¹³ as below:

${\dot{A}}^{\text{iARD}}=\dot{\gamma }[2D_r-2\text{tr}(D_r)A-5D_r\cdot A-5A\cdot D_r+10A_4:D_r]$

19

$D_r=C_I\left(I-C_M\frac{D^2}{\Vert D^2\Vert }\right)$

20

where D is the symmetric part of the velocity gradient tensor L, $D=\frac{1}{2}(L^{\text{T}}+L)$ . The scalar $\Vert D^2\Vert =\sqrt{\frac{1}{2}{D}^2:D^2}$ is the norm of tensor D².

Eventually, the RPR model is introduced as:

${\dot{A}}^{\text{RPR}}=-R\cdot {\dot{\Lambda }}^{\text{IOK}}\cdot R^{\text{T}}$

21

${\dot{\Lambda }}_{ii}^{\text{IOK}}=\alpha {\dot{\lambda }}_i,i,j,k=1,2,3$

22

where ${\dot{\Lambda }}^{\text{IOK}}$ is the material derivative of a particular diagonal tensor and its superscript indicates the intrinsic orientation kinetics assumption ^11–12 ; R is the rotation matrix and R^T is the transpose of R; the superscript T is the transpose operator of a matrix throughout this article; λ_i is the eigenvalues of A, ${\lambda }_1\ge {\lambda }_2\ge {\lambda }_3$ ; and $R=[e_1,e_2,e_3]$ is defined by eigenvector columns of A.

The iARD-RPR model has only the three physical parameters: a fiber–fiber interaction parameter C_I , a fiber–matrix interaction parameter C_M , and a slow-down parameter α. The available region of these three parameters is suggested in: 0 < C_I < 0.1, 0 < C_M < 1, 0 < α < 1. For short fibers, the default parameters are suggested: C_I = 0.005, C_M = 0.1, α = 0.7; for long fibers, C_I = 0.01, C_M =1.0, α = 0.1. In addition, the shell-orientation layer depends on the iARD parameters, while the core width is controlled by the RPR parameter. The value of fiber orientation in the shell increases with increasing C_I and decreasing C_M . The core width increases with increasing α. According to these criteria, we could fine-tune three parameters to obtain an accurate orientation prediction compared with experimental data, namely, determining an optimal parameters set.

Regarding a particular coupled problem, the flow and orientation equations must be solved simultaneously, since fibers orient in response to flow, and flow and suspension rheology are affected by fiber orientation in reality. To date, the coupled approach has not been implemented in commercial software packages due to require significant computational resources. Previously, VerWeyst and Tucker ³⁰ demonstrated that the effects of coupling the flow filed and fiber orientation calculations had a negligible effect on the solution. For engineering applications, a decoupled approach has been availably applied in which the fluid flow problem is first solved, and the resulting velocity field is used to compute the fiber orientation.

Results and discussion

Referring to the previous disclosed research report by the Pacific Northwest National Laboratory and the Oak Ridge National Laboratory ³¹ in the United States, as well as a previous study, ¹² a center-gated disk geometry is shown in Figure 1. The disk is about 90 mm in radius with a gap thickness of 2 h = 3 mm. One measurement point is located at 34 mm in the radial direction. The mold and melting temperatures were set at 70°C and 238°C, respectively. The filling time and volume flow rate were 0.59 s and 130 cc/s, respectively. The packing time of 10 s and packing pressure of 100 MPa were given, and all of the processing conditions are addressed in Table 1.

Figure 1.

The injection-molding geometry of the center-gated disk (90 mm radius and 3 mm thickness) and its hybrid mesh (20 layers) for 3D numerical computation.

Table 1.

Processing conditions for the 40 wt% LGF/PP composite during the injection mold filling.

Fill time	0.59
Fill rate (cc/s)	130
Mold temperature (°C)	85
Melt temperature (°C)	270
Pack time	10
Pack pressure	50

LGF/PP: long glass fiber immersed in polypropylene matrix.

The researched material is the 40 wt% long glass fiber immersed in polypropylene matrix (40 wt% LGF/PP) with a density of 1.22 g/cm³. The initial fiber diameter and length are 17.4 μm and 13 mm, respectively. The post-extruded average fiber length is about 2 mm.

All of material parameters, including the PVT, shear viscosity, specific heat, and thermal conductivity, are summarized in Tables 2 to 4.

Table 2.

The PVT model parameters of the 40 wt% LGF/PP composite.

b 1S	cm3/g	0.9781
b 1L	cm3/g·K	0.0005943
b 2S	dyne/cm2	1.2522 × 109
b 2S	1/K	0.005063
b 2L	cm3/g	0.9089
b 3S	cm3/g·K	0.0003226
b 3L	dyne/cm2	2.3069 × 109
b 4S	1/K	0.006511
b 5	K	425.67
b 6	cm2·K/dyne	6.6667 × 10−9
b 7	cm3/g	0.04196
b 8	1/K	0.1
b 9	cm2/dyne	6.6667 × 10−10

LGF/PP: long glass fiber immersed in polypropylene matrix.

Table 3.

The WLF–Cross viscosity model parameters of the 40 wt% LGF/PP composite.

N		0.2768
τ*	dyne/cm2	35,922
D 1	g/cm·s	3.3648 × 1017
D 2	K	263.15
D 3	cm2·K/dyne	0
A 1		38.39
∼A2	K	51.6
τ Y0	dyne/cm2	0.07
TY	K	9000

WLF: William–Landel–Ferry; LGF/PP: long glass fiber immersed in polypropylene matrix.

Table 4.

The specific heat and thermal conductivity constants of the 40 wt% LGF/PP composite.

Specific heat Cp (erg/g·°C)	2.64 × 107
K (erg/s·cm·°C)	17,300

LGF/PP: long glass fiber immersed in polypropylene matrix.

The 3D model was first developed by one of the Computer Aided Design programs and then meshed by creating hexagonal elements with 20 layers, as shown in Figure 1. In this study, the Moldex3D simulation program was used in the disk during the injection molding process. Figure 2 shows the sprue pressure history from the filling stage (0.59 s) to the packing stage (10 s).

Figure 2.

Sprue pressure history from filling stage to packing stage.

For fiber orientation prediction, the iARD-RPR model parameters are listed in Table 5. Additionally, the initial condition of the fiber orientation tensor at the runner entrance is set as the isotropic state. Due to the higher 3D numerical resolution and the runner system consideration in the computation, it was unnecessary to set an experimental inlet condition around a gate for enhancing the accuracy in determining the fiber orientation.

Table 5.

Optimal parameters of the iARD-RPR fiber orientation model for the 40 wt% LGF/PP composite.

CI	0.01
CM	1.0
α	0.1

iARD-RPR: Improved Anisotropic Rotary Diffusion and Retarding Principal Rate; LGF/PP: long glass fiber immersed in polypropylene matrix.

Both fiber orientation components, A ₁₁ and A ₂₂, correspond to the flow and the cross-flow directions, respectively. As a result, Figures 3 and 4 show the predicted gap-wise distribution of fiber orientation tensor components, A ₁₁ and A ₂₂, at the measurement point located at 34 mm in the radial direction. The solid and dash lines are at the end of the filling and packing stages, respectively. For the experimental data on fiber orientation, reference can be made to the PNNL–ORNL report ³¹ ; these predicted distributions exhibit a distinct core–shell structure. The fibers found in the shell region near the cavity are strongly aligned with the flow direction, as opposed to the others found in the core region near the cavity center, which are transverse to the flow. In particular, the core width of the orientation distribution for the filling is obviously broader than that for the packing. Significantly, this prediction is in qualitative agreement with the earlier experimental results on long fibers in the study by Bailey and Rzepka, ⁵ who concluded that the core thickness is reduced during packing. As a whole, the iARD-RPR predictions for the core region and shell layers agree fairly well with the PNNL–ORNL results of the fiber orientation data.

Figure 3.

Comparison of the predicted orientation distributions along the flow directions through normalized thickness at the measure point for the filling and packing stages with experimental validation for the 40 wt% LGF/PP composite. LGF/PP: long glass fiber immersed in polypropylene matrix.

Figure 4.

Comparison of the predicted orientation distributions along the cross-flow directions through normalized thickness for the filling and packing stages with experimental validation for the 40 wt% LGF/PP composite. LGF/PP: long glass fiber immersed in polypropylene matrix.

After a cavity has been completely filled, the density and viscosity of a compressible polymeric melt affect the state of fiber orientation during the packing stage. Bailey and Rzepka ⁵ proposed the mechanism of the large viscosity gradient from the lower viscosity of the core region to the higher viscosity of the shell layer during the packing stage. According to the simulation, Figure 5 shows the density profile at the measurement point through normalized thickness at this point, for a comparison of the filling and the packing stages. There is an obvious compressible effect during the packing stage since the packing density is higher than the filling density.

Figure 5.

Density profile through normalized thickness at the measure point for the filling and packing stages.

Moreover, Figure 6 presents the shear viscosity profile on a logarithmic scale for a comparison of the filling and packing stages. During the filling, the shell viscosity is lower than the core viscosity. In other words, fibers at the shell layer easily align along the flow direction due the lower flow resistance. Near the boundary close to the cold mold, a very high viscosity due to the frozen fluid is found. In contrast to the packing, the core viscosity is lower than the shell viscosity, namely involving an increase in the flow-direction fiber orientation A ₁₁. Thus, this predicted finding supports the viscosity gradient proposed by Bailey and Rzepka ⁵ for a reduction of core thickness or an increase of fiber orientation during the packing. This result is attributed to the packing-induced flow in the core region.

Figure 6.

Viscosity profile on a logarithmic scale through normalized thickness at the measure point for a comparison of filling stage and packing stage.

Furthermore, Figure 7 exhibits the flow-direction fiber orientation A ₁₁ distribution at two various packing times, 1.0 and 10 s under constant packing pressure of 100 MPa. Obviously, the A ₁₁ distribution is kept almost constant. In addition, there is a small change in the fiber orientation distribution at two packing pressures, 100 and 500 MPa under constant packing time of 10 s, as shown in Figure 8. There is a slight reduction of core thickness under a strong packing pressure. These predicted results also agree with Bailey and Rzepka’s experimental observation. ⁵

Figure 7.

Flow-direction fiber orientation composite normalized thickness at the measure point for various packing times.

Figure 8.

Flow-direction fiber orientation composite normalized thickness at the measure point for various packing pressures.

Conclusions

In the past, fiber orientation states were merely determined in the filling stage for injection molding simulation based on the 2.5D or Shell computational algorithm. In this study, we focused on the injection molding process of the center-gated disk from the filling to the packing stage for long glass fiber composite. By employing the state-of-the-art technology of 3D injection molding simulation, our ultimate goal was to demonstrate the effect of the packing stage on fiber orientation. The packing core width of the orientation distribution is obviously narrower than that of the filling due to the high viscosity gradient. Furthermore, the fiber orientation distribution is almost unchanged at various packing times, whereas the great packing pressure slightly influences the core thickness reduction. The predicted results also agree with the previous experiments. Therefore, we believe that the effect of the packing on fiber orientation is primary for a real, bigger, complex geometry of designing high-quality automotive parts, such as various changes in the direction of the flow, the inclusion of ribs, and the changes in the thickness and holes.