# Flow-fiber coupled viscosity in injection molding simulations of short fiber reinforced thermoplastics

Tianyi Li, Jean-François Luyé

► To cite this version:

Tianyi Li, Jean-François Luyé. Flow-fiber coupled viscosity in injection molding simulations of short fiber reinforced thermoplastics. International Polymer Processing, 2019, 34 (2), pp.158-171. hal-01683052v5

HAL Id: hal-01683052

https://hal.science/hal-01683052v5

Submitted on 15 Sep 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
T. Li*, J.-F. Luyé

Promold, Paris, France

# Flow-Fiber Coupled Viscosity in Injection Molding Simulations of Short Fiber Reinforced Thermoplastics

The main objective of this paper is to numerically investigate the use of fiber-dependent viscosity models in injection molding simulations of short fiber reinforced thermoplastics with latest commercial software. We propose to use the homogenization-based anisotropic rheological model to take into account flow-fiber coupling effects. The 4th-order viscosity tensor is approximated by an optimal scalar model and then implemented in the Moldflow Insight API framework. Numerical simulations are performed for a test-case rectangular plate with three fiber orientation models. The resulting coupled flow kinematics and fiber evolutions are then compared to the standard uncoupled simulations. Interpretations are given based on detailed post-processing of the field results. Certain deformation conditions are expected to be better taken into account, which may also in return lead to an improved fiber orientation prediction. Preliminary confrontation between flow-fiber coupled simulations and existing experimental data is then presented at the end of the paper.

## 1 Introduction

Short fiber reinforced thermoplastics are gaining popularity in industries because they can guarantee mechanical resistance specifications while achieving overall part weight reductions. Integrative simulation procedures (Adam and Assaker, 2014) are extremely appealing since heterogeneous local microstructure information caused by different processing steps can then be transferred to the subsequent mechanical simulations. In particular, fiber-induced material anisotropy wisely placed in some critical loaded regions can be better taken into account, resulting in an even optimized part design.

An experimentally validated injection molding simulation constitutes hence the cornerstone of these integrated simulations of injected short fiber reinforced thermoplastics. It is now well known that an inaccurate prediction of fiber orientation distribution in the part could affect ultimate structural response simulations via finite element analyses (Wedgewood

et al., 2017). According to Arif et al. (2014) and Rolland et al. (2016), a well predicted injection molding induced skin-shell-core fiber orientation structure ((Papathanasiou, 1997) for a comprehensive review on this subject) is also crucial since it contributes to a correct prediction of damaging mechanisms and ultimate failure of these composites.

It has been discovered that the fiber orientation models currently available in commercial software produce often unsatisfactory predictions in the core region (Kleindel et al., 2015; Tseng et al., 2017a). Compared to experimental results, fibers are estimated to be over-aligned in the flow direction on the mid-surface and the width of the core is also under predicted. Since the core region is in general dominated by extensional flows transverse to the filling direction, it can be argued that these elongational deformations are not correctly taken into account in current models and a flow-fiber coupled simulation with fiber-dependent rheological properties may be more appropriate (Jørgensen et al., 2017; Tucker III, 1991).

There exists already a relatively abundant literature on different flow-induced rheological models ((Phan-Thien and Zheng, 1997) for instance for a quick review). The basic idea is that under shear flows fibers contribute less to the overall suspension viscosity when they are aligned in the shear direction, while it is the opposite for elongational flows in the same direction (Laun, 1984). Since then these fiber-dependent viscosity models have been applied to flow and injection molding simulations (Chung and Kwon, 1996; Costa et al., 2015; Mazahir et al., 2013; Ranganathan and Advani, 1993; Redjeb et al., 2005; Tseng et al., 2017a; Verweyst and Tucker III, 2002; Vincent et al., 2005). The research is still on-going and no definite conclusions have been yet reached on whether or not these flow-fiber coupling effects are really important for most cases. That's probably one of the reasons why currently almost all commercial injection molding software still perform an uncoupled analysis.

The objective of this paper is thus to present some updated results on the subject with a latest injection molding software (Autodesk Moldflow Insight 2018, Autodesk Australia Pty Ltd., Kilsyth, VIC, Australia). Thanks to its API framework (Costa et al., 2015) in section 2 a representative anisotropic fiber-induced viscosity model is implemented as a user-defined viscosity function. In particular, fiber orientation is explicitly taken into account in the viscosity. Compared to some previous researches, in section 3 the flow-fiber coupling effects are numerically eval-

<footnote>* Mail address: Tianyi Li, Promold, 42 rue Boursault, 75017 Paris, France
E-mail: tianyikillua@gmail.com</footnote>
uated by using more physical injection conditions via non-isothermal three-dimensional finite element computations. Furthermore, besides the traditional Folgar-Tucker model (Folgar and Tucker III, 1984), two latest experimentally validated fiber orientation models (RSC and MRD) are also compared in terms of their possible different contributions to the overall coupling. A preliminary confrontation between the simulation results and existing experimental data available in the literature is then presented in section 4. Conclusions and future research directions drawn from the numerical results are indicated in section 5.

General notation conventions adopted in this paper are summarized as follows. Scalar-valued quantities will be denoted by italic Roman or Greek letters like temperature T and pressure p. Vectors and second-order tensors as well as their matrix representation will be represented by boldface letters. This concerns for example the 2nd-order fiber orientation tensor **a**, the stress tensor **σ** and the rate of deformation tensor **D**. Higher order tensors will be indicated by blackboard letters such as the 4th-order orientation tensor **A** and the 4th-order viscosity tensor **V**. Tensors are considered as linear operators and intrinsic notation is adopted. If the resulting quantity is not a scalar, the contraction operation will be written without dots, such as $(\nabla\mathbf{D})_{ij} = \nabla_{ijkl}\mathbf{D}_{kl}$. Inner products (full contraction giving a scalar result) between two tensors of the same order will be denoted with a dot, such as $\mathbf{D} \cdot \nabla\mathbf{D} = \nabla_{ijkl}\mathbf{D}_{kl}\mathbf{D}_{ij}$. The Einstein summation convention is assumed.

## 2 Flow-Fiber Coupled Viscosity Models

### 2.1 Anisotropic Fiber-Induced Rheological Equations

For (short) fiber reinforced thermoplastics (that is, when the immersed particle aspect ratio r is sufficiently large, say r > 20, which is generally satisfied by fibers), the anisotropy induced by these fibers on the overall suspension rheological properties can be described by the following expression of the stress tensor

$$ \sigma = -p\mathbf{I} + 2\eta\mathbf{D} + 2\eta N_p A \mathbf{D}, \quad (1) $$

where **σ** is the (macroscopic) stress tensor for the fiber reinforced suspension, **p** is the pressure, **I** is the 2nd order identity tensor, **D** = $\frac{1}{2}\text{dev}(\nabla\mathbf{v} + \nabla^T\mathbf{v})$ is the deviatoric deformation rate tensor, $\eta$ refers to the viscosity of the matrix (without fibers) and **A** designates the 4th order fiber orientation tensor introduced in (Advani and Tucker III, 1987). The interested reader can find detailed theoretical explanation of Eq. 1 and homogenization-based derivation of similar models in Lipscomb et al., 1988; Phan-Thien and Zheng, 1997; Tucker III, 1991) and references therein. Note that according to Eq. 1, the viscosity of the suspension is now no longer a scalar but a 4th order tensor reflecting the anisotropy induced by fibers

$$ \sigma = -p\mathbf{I} + 2\nabla\mathbf{D} \quad \text{with} \quad \nabla = \eta(\mathbf{I} + N_p A), \quad (2) $$

with **I** the 4th order identity operator. In Eq. 1 and 2, the coefficient $N_p$ measures the scalar intensity of the anisotropic contribution of fibers to the overall viscosity. In particular when one sets $N_p = 0$, we retrieve the uncoupled scalar viscosity model, completely independent of fiber evolutions during the injection process. Since traditional linear dependence of $N_p$ on volume

fraction may not be suitable for concentrated fiber-reinforced suspensions (Jørgensen et al., 2017), in this paper we adopt the nonlinear functional dependence of $N_p$ on fiber aspect ratio and volume fraction proposed by Phan-Thien and Graham (1991). Its dependence on mass fraction is illustrated in Fig. 1 for a typical 50%wt fiber filled polyamide (PA-GF50) with the fitting parameter $A = 50\%$. For comparison, the dilute theory of Ericksen (1959) is also indicated. These two theories agree for small fiber volume concentrations, however for the current case (50%wt), the nonlinear model of Phan-Thien and Graham (1991) predicts a much larger anisotropic contribution coefficient $N_p$.

### 2.2 Optimal Scalar Approximation of the 4th-Order Viscosity Tensor

In order to evaluate numerically the anisotropic fiber-induced viscosity model in injection molding simulations, we propose to implement Eq. 2 in the Moldflow Insight API framework described in Costa et al. (2015).

Currently, only a scalar user-defined viscosity function is supported via this interface. Hence, the original model as shown in Eq. 2 needs to be adapted for implementation. One novelty of this paper consists in proposing an optimal approximate scalar (isotropic) viscosity value of the anisotropic 4th order model.

According to Eq. 2, for stress computations the viscosity tensor **V** is used via its action $\nabla\mathbf{D}$ on the current deformation rate tensor **D**. The idea is thus to define an effective scalar viscosity value $\eta_*$ such that the scalar multiplication $\eta_*\mathbf{D}$ is as close as possible to $\nabla\mathbf{D}$ in a certain sense. In this paper, we simply use the classical Frobenius norm (square root of the sum of its components squared) for 2nd order tensors and we minimize the following approximation error

$$ \|\eta_* \mathbf{D} - \nabla\mathbf{D}\|^2 = \min_{\mathbf{v}} \|\mathbf{v}\mathbf{D} - \nabla\mathbf{D}\|^2. \quad (3) $$

Through the Frobenius (elementwise) inner product, the scalar $\eta_*$ can be simply regarded as the projection of **V** in the direction of **D**. By a direct calculation of Eq. 3, we obtain thus the

<figure><img src="image_2.png" /><figcaption>Fig. 1. Anisotropic contribution coefficient $N_p$ for PA-GF50 as a function of the mass fraction using the dilute (Ericksen, 1959) and the concentrated (Phan-Thien and Graham, 1991) theories</figcaption></figure>
optimal scalar viscosity approximating the 4th-order anisotropic viscosity when $D \neq 0$

$$ \eta_{*} = (1 + N_{p}a_{*})\eta \quad \text{with} \quad a_{*} = \frac{D \cdot A D}{\|D\|^{2}} = \frac{\text{tr}(D^{\text{T}}A D)}{\|D\|^{2}}, \quad (4) $$

where tr is the trace operator for 2nd-order tensors. In Eq. 4, the scalar $a_*$ can also be regarded as the optimal equivalent scalar of the 4th order orientation tensor A, for the current deformation rate tensor D considered. For $D = 0$, we can simply set $a_* = 0$.

The computation of this optimal scalar viscosity through Eq. 4 requires the 4th order fiber orientation tensor A. In general in injection molding simulations it can only be recovered approximately from the 2nd order tensor via a particular closure formulation. For consistency, we use the same orthotropic (ORT) closure model described in Verweyst (1998) which is used by default in Moldflow since the version 2017R2. The tensor A is positive semi-definite by definition (Advani and Tucker III, 1987). Accordingly the scalar $a_*$ is non-negative and the optimal scalar suspension viscosity $\eta_*$ remains equal to or greater than the matrix viscosity $\eta$. It is desired that the closure formulation also satisfies this property. In Fig. 2, the smallest eigenvalue of A approximated by the ORT closure is numerically computed for all possible 2nd-order tensor orientations with $(a_1, a_2)$ the two largest eigenvalues of a, in the TUB triangular domain of Cintra Jr. and Tucker III (1995). It can be seen that the ORT-approximated A only loses definiteness for planar or unidirectional orientation states (generate cases of 3d orientation states). The semi-positiveness of the ORT model is thus numerically verified.

In order to better understand the optimal scalar viscosity model Eq. 4, its behavior under two typical flows (simple shear and planar elongation) frequently present in injection molding is considered and then compared to the original model (Eq. 2). The rate-of-deformation tensor D corresponding to these two flows is respectively given by

$$ D_1 = \frac{1}{2} \begin{bmatrix} 0 & 0 & \dot{\gamma} \\ 0 & 0 & 0 \\ \dot{\gamma} & 0 & 0 \end{bmatrix}, \quad D_2 = \begin{bmatrix} 1 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & -1 \end{bmatrix} \dot{\epsilon}. \quad (5) $$

<figure><img src="image_3.png" /><figcaption>Fig. 2. Smallest eigenvalue of the 4th-order fiber orientation tensor approximated by the orthotropic closure formulation in the TUB orientation space</figcaption></figure>

In the case of a simple shear in the 1–3 plane given by $D_1$ in Eq. 5, a direct application of Eq. 4 along with the stress expression (Eq. 2) by replacing the 4th-order viscosity V by its scalar approximation $\eta_*$ gives

$$ \begin{align} a_* &= 2A_{1313} \Rightarrow \sigma_* \nonumber \\ &= \begin{bmatrix} -p & 0 & (1+2N_pA_{1313})\eta\dot{\gamma} \\ 0 & -p & 0 \\ (1+2N_pA_{1313})\eta\dot{\gamma} & 0 & -p \end{bmatrix}. \quad (6) \end{align} $$

It can be seen that the most adequate scalar value of the 4th-order fiber orientation tensor for simple shear flows according to Eq. 4 is its component in the shear plane. Comparing the approximated stress tensor (Eq. 6) to the original one given by Eq. 1, one finds that the 13 stress component is exactly the same. The shear viscosity is thus given by

$$ \eta_s = \sigma_{13}/\dot{\gamma} = (1 + 2N_p A_{1313}) \eta, \quad (7) $$

for both the original and optimal scalar models, and however as expected the approximate model does not present any normal stress differences, while for the original one we have for instance $N_1 = \sigma_{11} - \sigma_{22} = 2N_p(A_{1113} - A_{2213})\eta\dot{\gamma}$.

As for planar extensional flows given by $D_2$ in Eq. 5, it can also be easily shown that the planar extensional viscosity (Petrie, 2006) is exactly recovered during the scalar approximation process (Eq. 4), which is given by

$$ \begin{align} \eta_p &= \frac{\sigma_{11} - \sigma_{33}}{\dot{\epsilon}} = 2(2 + N_p A_{1111} - 2N_p A_{1133} + N_p A_{3333})\eta \nonumber \\ &= 4(1 + N_p a_*)\eta. \quad (8) \end{align} $$

One can check that when fiber contributions ($N_p$ or $A_{ijkl}$) are small, one retrieves the theoretical value $\eta_p = 4\eta$ for a Newtonian fluid (Petrie, 2006).

As an illustration of the fiber orientation effects present in the optimal scalar model (Eq. 4), we consider two particular orientation states: a fully random state $a_R$ and a quasi-perfect unidirectional orientation $a_{UD}$. The 2nd-order fiber orientation tensors are respectively given by

$$ \begin{align} a_R &= \frac{1}{3} \mathbf{I} = \frac{1}{3} \begin{bmatrix} 1 & & \\ & 1 & \\ & & 1 \end{bmatrix}, \\ a_{UD} &= \begin{bmatrix} \bar{a}_{11} & & \\ & \frac{1}{2}(1 - \bar{a}_{11}) & \\ & & \frac{1}{2}(1 - \bar{a}_{11}) \end{bmatrix}, \quad (9) \end{align} $$

with $\bar{a}_{11} = 0.98$. In Fig. 3, the values of the coupling factor $a_*$ for these two particular orientation states under two typical flows (Eq. 5) are indicated. We can see that under simple shear, the coupling factor, and hence also the shear viscosity according to Eq. 6, are larger for the random orientation state than the unidirectional one. In the latter case the fiber contribution is negligible. For the planar extensional case, fibers contribute the most to the planar extensional viscosity (Eq. 8) when they are aligned with the extension direction. Their contribution is reduced when they are randomly oriented in space. These theoretical observations are well conforming to the experimental findings presented for instance in Laun (1984).
<figure><img src="image_4.png" /><figcaption>Fig. 3. Coupling factor $a_*$ for two particular orientation states under two typical flows, for UD: $\bar{a}_{11} = 0.98$</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 4. Variation of $\bar{a}_*$ as $\bar{a}_{11}$ increases from 1/3 (isotropic state) to 1 (unidirectional state)</figcaption></figure>

## 2.3 Implementation in the Moldflow Insight API Framework

We recall that in Eq. 4, the viscosity $\eta$ refers to that of the matrix in absence of fibers. In Moldflow API simulations, it is more convenient to still select the fiber reinforced material (say 50%wt) with a measured viscosity $\eta_{MF}$, based on which we will then apply a scaling factor $\hat{\eta}$ incorporating anisotropic effects according to Eq. 4. Consequently, we need to recover from $\eta_{MF}$ at least approximately the viscosity function $\eta$ corresponding to the non-filled matrix. Since in Moldflow the apparent viscosity $\eta_{MF}$ is measured using a capillary rheometer at relatively high shear rates, the fibers should be almost perfectly aligned in the shear direction (Costa et al., 2015). We consider hence a simple shear flow $v = x_3 \dot{\gamma} e_1$ in the 1-3 plane where the deformation rate is given by $D_1$ in Eq. 5. Due to possible fiber-fiber/matrix interactions (Folgar and Tucker III, 1984, for instance), the alignment is not perfect and we suppose that the 2nd order fiber orientation tensor can be parameterized by its 11 component and is given by $a_{UD}$ in Eq. 9, where $\bar{a}_{11}$ should be in general very close to 1 characterizing a quasi unidirectional state. Applying our optimal scalar viscosity (Eq. 4) for $a_{UD}$ under the simple shear flow gives

$$ \eta_{MF} = (1 + N_p \bar{a}_*) \eta, \quad (10) $$

where $\bar{a}_*$ is obtained by combining Eq. 4 with $a_{UD}$ in Eq. 9 and $D_1$ in Eq. 5

$$ a_* = \frac{D_1 \cdot A D_1}{\|D_1\|^2}, \quad A = A(a_{UD}(\bar{a}_{11})). $$

Its functional dependence with respect to $\bar{a}_{11}$ is illustrated in Fig. 4. We verify that for $\bar{a}_{11} = 1$ when the fibers are (somehow indeed) perfectly aligned in the shear direction, the Moldflow apparent viscosity $\eta_{MF}$ directly gives the matrix viscosity $\eta$ since in this case fibers have no effect on the suspension viscosity, resulting in $\bar{a}_*$ = 0. For other cases, this small coefficient $\bar{a}_*$ reflects thus the fiber alignment information contained in the Moldflow fiber reinforced suspension viscosity. Two particular values of $\bar{a}_{11}$ that will be used in the following numerical simulations in section 3 are also indicated in Fig. 4.

Combining Eq. 4 and Eq. 10, we finally obtain the optimal scalar fiber-induced rheological model that will be implemented in the Moldflow API framework

$$ \begin{aligned} \eta_*(\mathbf{D}, \mathrm{T}, \mathrm{p}, \mathbf{a}) &= \hat{\eta}(\mathbf{D}, \mathbf{a}) \cdot \eta_{\mathrm{MF}}(\dot{\gamma}, \mathrm{T}, \mathrm{p}) \\ &= \frac{1 + N_p \bar{a}_*(\mathbf{D}, \mathbf{a})}{1 + N_p \bar{a}_*} \eta_{\mathrm{MF}}(\dot{\gamma}, \mathrm{T}, \mathrm{p}). \end{aligned} \quad (11) $$

We recall that here $\hat{\eta}$ is the non-dimensional viscosity scaling factor that will be applied to the Moldflow fiber reinforced suspension viscosity $\eta_{MF}$. Observe that now the shear rate (tensor) dependence of the resulting semi-anisotropic fiber dependent viscosity $\eta_*$ is both through the original Moldflow viscosity function (mainly the Cross model) and the coefficient $a_*$ adapting fiber orientation $a$ induced anisotropy to the current deformation rate tensor $D$ via Eq. 4. The temperature $T$ and pressure $p$ dependence remain unchanged compared to standard un-coupled simulations.

Finally, from a purely algorithmic point of view, the above orientation-dependent rheological equation is implemented under Moldflow Insight using a weakly-coupled approach according to Laure et al. (2011) and Redjeb et al. (2005). This means at each time-step $t^n$ the flow and the fiber equations are solved sequentially and independently, and the fiber orientation used to calculate the current viscosity $\eta_*(D^n, T^n, p^n, a^{n-1})$ in Eq. 11 is taken from the last time-step $t^{n-1}$. Admittedly, a strongly-coupled numerical implementation (Ranganathan and Advani, 1993; Verweyst and Tucker III, 2002)) should be employed to ensure that the flow and the fiber solutions are always in phase. Unfortunately it is not (yet) possible via the current Moldflow Insight API framework. Nevertheless, it can be expected that by reducing the time-step increment, the “phase-shift” between flow and fiber will be reduced, converging thus to a strongly coupled solution. For readability, a qualitative illustration of this point based on the numerical results to be presented in section 3 is given here. Under Moldflow, the time-step increment can be controlled by the easily understandable “maximum %volume to fill per time step” parameter $\Delta V$. Its influence on fiber orientation in our weakly-coupled setting is presented in Fig. 5 for three different $\Delta V$s. No numerical divergence issues are reported during all of these simulations.
<figure><img src="image_6.png" /><figcaption>Fig. 5. Fiber orientation $a_{22}$ on the midsurface along the length direction from the gate $x_1/L = 0$ up to the end $x_1/L = 1$ obtained for three different time-step increments, $z = 0$</figcaption></figure>

tions. Surprisingly, the fiber orientation results remain quite insensitive to this parameter except near the end of fill. This convergence behavior should depend on the exact process settings and mostly on the particular material properties. In this case it can thus be considered that the default $\Delta V = 4\%$ provides already a quite good approximation to the strongly coupled solution. Note that the subsequent numerical simulations are conducted with $\Delta V = 1\%$.

### 3 Test-Case Simulations on a Rectangular Plate

In this section we intend to investigate numerically the effect of fiber-induced viscosity models presented in section 2 through injection molding simulations under Autodesk Moldflow. The version Moldflow Insight 2018 is used. The Moldflow study (.sdy) file is kindly provided by the first author of Costa et al. (2015) and we use the same geometry, basic material, injection process and numerical parameters as used in Costa et al. (2015). As can be seen from Fig. 6A, the geometry consists of a simple end-gated rectangular flat plate with a thickness of 3.6 mm. Two probe points P1 (upstream, near the gate) and P2 (downstream, at the plate center) are used to analyze the simulation. The default material parameters of a 50%wt short fiber filled polyamide (Zytel PLS95G50DH2 BK261) in the Moldflow database are used, with a standard Cross-WLF viscosity model. The $N_p$ parameter of Phan-Thien and Graham (1991) is computed using the aspect ratio and volume fraction information contained in the Moldflow database, and a fitting parameter $A = 50\%$, which gives $N_p \approx 221$ (see also Fig. 1). Fiber orientation may continue to evolve during the subsequent packing phase, however for simplicity the present investigation is limited to the filling phase. For comparison concerns we fix a constant ram speed at 50 mm/s which leads to an injection time near 1.15 s. Spatial discretization is realized with sufficiently fine triangles on the surfaces and 16 layers of tetrahedral elements in the gapwise direction.

It will be interesting to evaluate the flow-fiber coupling effects introduced by the fiber-induced viscosity models presented in section 2, for several fiber orientation models which may contribute differently to the overall coupling. In this paper

<figure><img src="image_7.png" /><figcaption>Fig. 6. Rectangular end-gated flat plate with two probe points (A), filling contour for the uncoupled and coupled injection simulations (B)</figcaption></figure>

we will use the three principal short fiber orientation models available in Moldflow Insight 2018: the classical Folgar-Tucker model (Folgar and Tucker III, 1984), the Reduced Strain Closure (RSC) model proposed by Wang et al. (2008) and the lately introduced Moldflow Rotational Diffusion (MRD) model (Autodesk, 2016). Note that since the Moldflow version 2017R2 the default fiber orientation model becomes the MRD model. The default parameters for these fiber models provided by Moldflow are used throughout the paper, which gives respectively

* Folgar-Tucker: $C_i = 0.0102$.

* RSC: $C_i = 0.002$ and $\kappa = 0.1$.

* MRD: $C_i = 0.0015$, $D_1 = 1$, $D_2 = 0.8$ and $D_3 = 0.15$.

#### 3.1 Fiber Orientations Predicted by Uncoupled Simulations

We will first compare these three fiber orientation models using the default uncoupled viscosity model (Cross-WLF). Since the flow equations and the fiber orientation model are uncoupled, the same filling pattern is obtained as illustrated in Fig. 6B.

The steady $a_{11}$ and $a_{13}$ fiber orientation components in the thickness direction of the probe point P1 are presented in Fig. 7 for these three models. At P2 similar results can be obtained and hence are not illustrated here. We recall that the 1-axis corresponds to the horizontal flow direction while the 3-axis is the thickness one, see Fig. 6A. The thickness coordinate is normalized such that $z = -1$ (resp. $z = +1$) corresponds to the lower (resp. upper) surface and $z = 0$ designates the mid-surface. As for the $a_{11}$ component, the well-known (skin)-shell-core structure can be observed for these three fiber mod-
<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 7. Fiber orientation $a_{11}$ (A) and $a_{13}$ (B) at P1 for the Folgar-Tucker, RSC and MRD models without flow-viscosity coupling</figcaption></figure>

els. The Folgar-Tucker model predicts a more flow-aligned orientation in the core as well as on the skin layer, whereas the other two models are designed to capture slower orientation kinetics observed in experiments. These two models also predict a reduced fiber orientation growth rate for the $a_{13}$ component especially in the shell and skin regions.

## 3.2 Coupled Simulations Using the Folgar-Tucker Model

The classical Folgar-Tucker model is first used to evaluate possible flow-fiber coupling effects introduced by Eq. 11. The $\bar{a}_{11}$ parameter characterizing fiber alignment in the Moldflow viscosity (Fig. 4) is equal to 0.93. The fiber orientation in the flow direction $a_{11}$ is presented in Fig. 8. It can be seen that while the orientation-dependent model indeed predicts a slightly decrease in the $a_{11}$ component at P1 especially in the core region, the coupling effect is not so obvious as can also be seen at P2. The flow-fiber coupled viscosity model may not be enough to predict a general slowed down fiber alignment, probably because fiber orientation predicted by the Folgar-Tucker model evolves too fast for any kind of flow kinematics. The flow-fiber coupling effects may become more visible and important when using a fiber model that explicitly and directly slows down general orientation kinetics observed in experiments.

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 8. Fiber orientation $a_{11}$ at P1 (A) and P2 (B) probe points using the Folgar-Tucker model with or without flow-fiber coupling where $\bar{a}_{11} = 0.93$</figcaption></figure>

## 3.3 Coupled Simulations Using the MRD Model

We then test the optimal scalar fiber-induced rheological model (Eq. 11) using the MRD fiber model, with as before $\bar{a}_{11} = 0.93$. The $a_{11}$ variation in the thickness for both probe points are presented in Fig. 9. At P1, the flow-fiber coupled model predicts a widened fiber core region compared to the uncoupled case due to the flow-fiber coupling. It also shows a slightly more oriented shell region. At P2, fibers are somehow predicted to be more aligned in the flow direction compared to the uncoupled model. Nevertheless the fiber orientation still displays a plausible skin-shell-core structure. With a slowed-down orientation growth rate, the MRD model displays more visible flow-fiber coupling effects compared to the previous Folgar-Tucker model. Note that a slower fiber alignment in the core region shown in Fig. 9A has also been observed with various viscosity models that take into account fiber orientation effects, cf. the references mentioned in section 1. On the other hand, depending on the circumstances, flow-fiber coupling may also lead to a slightly increased fiber alignment as observed here in Fig. 9B, cf. Laure et al. (2011).

The velocity profile and viscosity at P1 are then analyzed in Fig. 10 to evaluate the effect of flow-fiber coupling on kinematics and rheological properties. It can be seen from Fig. 10A
<figure><img src="image_12.png" /></figure>

A)

<figure><img src="image_13.png" /></figure>

B)

Fig. 9. Fiber orientation $a_{11}$ at P1 (A) and P2 (B) probe points using the MRD model with or without flow-fiber coupling where $\bar{a}_{11} = 0.93$

<figure><img src="image_14.png" /></figure>

A)

<figure><img src="image_15.png" /></figure>

B)

Fig. 10. Results obtained at the end of filling at P1 using the MRD model with or without flow-fiber coupling where $\bar{a}_{11} = 0.93$: (A) velocity profile $v_x$; (B) viscosity profile

that compared to the uncoupled case, the velocity profile is now essentially flattened (plug flow) in the core region resulting from the corresponding viscosity increase (by a factor of 10 at the center) as indicated in Fig. 10B. The overall deformation rate is thus significantly reduced in that region, which explains a less aligned core in Fig. 9A. A similar plug flow has also been reported by Tseng et al. (2017a), where a yield stress was added into the standard Cross viscosity model in order to account for fiber orientation effects at lower shear rates. Since the probe point P1 is near the gate, we suspect that this fiber-induced increased viscosity should originate from an elongational flow in the transverse direction. This point will be verified in the next section.

Concerning the overall filling pattern when using the flow-fiber coupled model (Eq. 11), it is essentially the same as that illustrated in Fig. 6B for the uncoupled case. Similar to the orientation-dependent rheological model proposed in Costa et al. (2015), the viscosity increase only concerns a limited core region (according to Fig. 10B, only about 25% of the thickness). That's likely to be the reason why it fails to predict a locally concave contour reported in Akay and Barkley (1993), where the flow advances faster near the edge than in the center.

## 3.4 Coupled Simulations with the RSC Model

Now the flow-fiber coupling effects introduced by Eq. 11 will be evaluated for the RSC model. This time, two values (0.93 and 0.86) of the $\bar{a}_{11}$ parameter are used. Again we first compare the fiber orientation profile at both probe points for the coupled and uncoupled cases. According to Fig. 11, fiber orientation is not very sensitive to the $\bar{a}_{11}$ parameter. Interestingly, compared to the MRD case (Fig. 9), now for both probe points we have a slowed down fiber alignment throughout the thickness and consequently a widened fiber core region. This indicates that different fiber models may also contribute differently to the overall coupling effect.

In Figs. 9 and 11, only the final fiber orientation states at the end of filling are presented. It will also be interesting to analyze temporal fiber evolutions for both coupled and uncoupled models. The averaged (in the thickness direction) fiber orientation as a function of filling time is illustrated in Fig. 12. It can be observed that the temporal fiber orientation increases monotonically. Due to flow-induced viscosity increase, the overall deformation rate is reduced implying a systematic slowed down fiber alignment during the entire filling process.
<figure><img src="image_16.png" /></figure>

A)

<figure><img src="image_17.png" /></figure>

B)

<caption>Fig. 11. Fiber orientation $\bar{a}_{11}$ at P1 (A) and P2 (B) using the RSC model with or without flow-fiber coupling</caption>

The spatial variation of the fiber orientation tensor at the end of filling can be used to get a bigger picture of the coupling effect. The $\bar{a}_{11}$ and $\bar{a}_{22}$ components along the centerline indicated in Fig. 6A from the injection gate to the end of the plate are presented in Fig. 13 for two height levels. The curve at $z = 0$ represents the flow-induced fiber alignment on the mid-surface (thus in the core region), while the curves at $z = +0.8$ and $z = -0.8$ are averaged into a single curve $z = \pm 0.8$ characterizing the spatial variation of fiber orientation in the shell region. The $x_1$ distance from the injection gate is normalized by the length L of the plate, and two probe points are also shown.

* In the shell region dominated by high shear flows, the fibers are quickly oriented parallel to the flow direction $e_1$ near the gate and then remain relatively aligned ($\bar{a}_{11}$ is between 0.7 and 0.8) throughout the length direction except for the last 10% portion just before the end of filling, where a sudden drop in $\bar{a}_{11}$ as well as an increase in $\bar{a}_{22}$ are observed due to boundary effects.

* In the core region, the fibers are initially more oriented in the transverse direction probably due to high extensional flows toward the $e_2$ direction near the gate. In the central part of the plate 20% $\le \frac{x_1}{L} \le 80\%$ while fibers are advected by the filling process, a gradual increase of the $\bar{a}_{11}$

<figure><img src="image_18.png" /></figure>

A)

<figure><img src="image_19.png" /></figure>

B)

<caption>Fig. 12. Temporal evolution of the average $\bar{a}_{11}$ in the thickness at P1 (a) and P2 (b) using the RSC model with or without flow-fiber coupling where $\bar{a}_{11} = 0.93$</caption>

value can be observed and the fibers become slightly more aligned with the flow direction.

* As expected (Fig. 3), the flow-fiber coupling effects are more visible for the core region where fibers are not highly aligned with the flow. The orientation-induced viscosity increase leads to an overall slowed down fiber alignment in that region, throughout the centerline direction. On the contrary, the use of an orientation-dependent viscosity is insignificant in the shell region, since fibers are mostly oriented parallel to the dominant shear direction. These observations agree with the previous study of Mazahir et al. (2013) on a center-gated disk.

Now, we will investigate what exactly happened at these two probe points P1 and P2 that leads to what we observe in the previous figures. For that, we need to analyze the deformation rate tensor $D$ in the thickness direction. Since the tensor components of $D$ are not directly available in Moldflow, the velocity field $v$ at the end of fill is exported through a home-made script to another post-processing software, where the gradient of $v$ can be automatically computed. To make sure that our computation of $D = \frac{1}{2} \text{dev}(\nabla v + \nabla^T v)$ matches that of Moldflow, the generalized shear rate $\dot{\gamma} = \sqrt{2}\|D\|$ is first compared in Fig. 14. The shear rate given by Moldflow can be indeed accurately re-
covered, which validates our post-processing approach. The small discrepancy mainly originates from different interpolation schemes of the velocity and averaging processes.

It can be seen from Fig. 15A that the 13 shear component is dominant compared to other shear components, which is fairly typical for thin-walled injection molding (Tucker III, 1991).

<figure><img src="image_20.png" /><figcaption>Fig. 13. Spatial variation of the $a_{11}$ (A) and $a_{22}$ (B) components along the centerline from the injection gate to the plate end for two height levels: $z = 0$ (mid-surface, core region) and $z = \pm 0.8$ (shell region). The coupled results are obtained using the RSC model and $\bar{a}_{11} = 0.93$. Solid line: uncoupled, dashed line: coupled</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>Fig. 14. Reconstruction of the generalized scalar shear rate at P1</figcaption></figure>

According to Fig. 15B, a planar elongational flow toward the transverse direction is indeed observed especially on the mid-surface $z = 0$ at P1 near the injection gate, where the deformation rate tensor is given by

$$ \mathbf{D} = \begin{bmatrix} -4.9 & 0.8 & 0.9 \\ 0.8 & 5.3 & 0 \\ 0.9 & 0 & -0.4 \end{bmatrix} \text{1/s.} $$

Using the reconstructed deformation rate tensor in Fig. 15 and the fiber orientation tensors in Fig. 11A, the analytical formula of the optimal scalar flow-fiber coupled model (Eq. 11) can be applied. The thus obtained theoretical viscosity scaling factor is compared to that produced by our Moldflow API user routine in Fig. 16A. A relatively good agreement can be found between them, which hence validates our implementation of (Eq. 11) in the Moldflow Insight API framework. Again, the slight discrepancy may mainly originate from the exact calculation of the deformation tensor and especially different averaging operators in Moldflow and our post-processing. On the mid-surface a viscosity scaling factor near 10 is retrieved (Fig. 10). This illustrates that it is indeed the planar extensional flow in the 2 direction together with a nearly transverse fiber orientation state in the core region as shown in Fig. 11 that are

<figure><img src="image_22.png" /><figcaption>Fig. 15. Deformation rate at P1 using the RSC model with flow-fiber coupling where $\bar{a}_{11} = 0.93$: (A) shear components and (B) elongational components</figcaption></figure>
<figure><img src="image_23.png" /><figcaption>A)</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>B)</figcaption></figure>

<caption>Fig. 16. Viscosity scaling factor with $\bar{a}_{11} = 0.93$: (A) comparison at P1 between theoretical and Moldflow API user routine values using the RSC model; (B) comparison between the MRD and RSC fiber models at P2 (API results)</caption>

mainly responsible for the viscosity increase in the core region, conforming to our previous theoretical analysis in Fig. 3.

The above analyses are now applied to the probe point P2 in Fig. 17. This time the deformation rate tensor on the mid-surface is now given by

$$ \mathbf{D} = \begin{bmatrix} -1.7 & 0.1 & -0.2 \\ 0.1 & 0.7 & 0.1 \\ -0.2 & 0.1 & 1.1 \end{bmatrix} \text{1/s.} $$

It can be observed that at the downstream point P2 we no longer have a planar divergent flow but a biaxial extensional flow toward the 2-3 directions in the core region. Furthermore, in terms of magnitude it is also not comparable to that observed at the upstream probe point P1. Besides the advection effect illustrated in Fig. 13, there could be something else that maintains a widened (less aligned) fiber orientation core observed with the RSC model at P2 and also throughout the length direction, see Fig. 11B and Fig. 13. In Fig. 16B, we compare the viscosity scaling factor obtained with the MRD and the RSC model. Although on the mid-surface the MRD model gives a slightly larger value, the RSC model presents a larger overall value in the whole thickness. Since the 13 shear component is dominant almost throughout the thickness except very near

<figure><img src="image_25.png" /><figcaption>A)</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>B)</figcaption></figure>

<caption>Fig. 17. Deformation rate at P2 using the RSC model with flow-fiber coupling where $\bar{a}_{11} = 0.93$: (A) shear components and (B) elongational components</caption>

the mid-surface, the flow-fiber coupling should mainly come from the fibers oriented in that plane. According to Fig. 7B, by default the RSC model predicts indeed more fiber orientations in the 13 plane than the MRD model. In the authors' opinion this maintains an overall viscosity comparable to the nominal value used in the uncoupled analysis. When combined with the advection effect under an upstream elongational flow, it then slows down fiber orientation in the core due to a viscosity scaling factor still larger than 1 in that region.

Finally, in Fig. 18, we compare the pressure temporal evolution at the injection point, for two values of $\bar{a}_{11}$. A larger value (closer to 1) of $\bar{a}_{11}$ leads to an overestimated evolution compared to the uncoupled case. This agrees well with our definition in (Eq. 10) since in the limiting case $\bar{a}_{11} = 1$, we have $\bar{a}_* = 0$ and the uncoupled Moldflow viscosity is assumed to be that of the unfilled matrix, based on which we will apply then a scaling factor $1 + N_{pa_*}$ always bigger than 1. This may very likely increase the overall pressure prediction. The use of a smaller $\bar{a}_{11}$ assumes that fibers are not perfectly aligned during the Moldflow measurement of viscosity, and compensates the pressure prediction. This parameter $\bar{a}_{11}$ can thus be regarded as a fitting parameter of the model.
<figure><img src="image_27.png" /><figcaption>Fig. 18. Pressure at the injection point using the RSC model with flow-fiber coupling</figcaption></figure>

## 4 Preliminary Confrontation with Experimental Data

In this section we propose a preliminary confrontation between flow-fiber coupled simulation results and the experimental data first presented in Tseng et al. (2018). For this, we consider the injection of a 50%wt short-glass-fiber reinforced PA66 (Ultra-mid A3WG10, BASF, Wyandotte, MI, USA) into a simple fan-gated rectangular plate with one measurement point B located at the center, see Fig. 19. The basic material parameters from the Moldflow database are used. For flow-fiber coupled simulations, the N<sub>p</sub> parameter of Phan-Thien and Graham (1991) is again computed with a fitting parameter A = 50%, which gives N<sub>p</sub> ≈ 163. The process conditions are recalled in Table 1. As in Tseng et al. (2017b; 2018), only the filling phase is considered and the default “fibers aligned at skin/random at the core” inlet condition is prescribed at the injection point. Spatial discretization is performed with sufficiently fine tetrahedral elements with 12 layers in the thickness direction. Simulations conducted on another 20 layers-based mesh do not present any particular differences.

For this preliminary experimental validation, the RSC fiber orientation model is used with its default parameters C<sub>i</sub> = 0.002 and κ = 0.1. The ā<sub>11</sub> parameter as shown in Fig. 4 is equal to 0.98 for coupled simulations. Note that this particular value of C<sub>i</sub> is probably not an optimized choice for the material considered. However in this paper the focus is placed on flow-fiber coupling effects. The interaction between different fiber orientation parameters and the proposed coupled viscosity model will be investigated in our future work.

The gapwise variations of the a<sub>11</sub> and a<sub>22</sub> components of the fiber orientation tensor at the central probe point B are presented in Fig. 20, where the experimental data used in Tseng et al. (2018) are also indicated.

It can be seen that the simulated fiber orientations present an obvious skin-shell-core structure as confirmed by the experiments. Furthermore, with the default parameters the simulation results agree quite well with the measured data in particular on the mid surface (z = 0) and for the a<sub>22</sub> component. Similar to Fig. 11, the use of a flow-fiber coupled viscosity slows down a little bit the general fiber alignment and decreases further the discrepancy between predicted and experimental results.

<figure><img src="image_28.png" /><figcaption>Fig. 19. Geometry of a fan-gated rectangular plate with one probe point located at the center</figcaption></figure>

<table><caption>Table 1. Process conditions for the used material</caption><thead><tr><td>Injection time s</td><td>Mold surface temperature °C</td><td>Melting temperature °C</td></tr></thead><tbody><tr><td>1.8</td><td>85</td><td>290</td></tr></tbody></table>

<figure><img src="image_29.png" /></figure>

A)

<figure><img src="image_30.png" /></figure>

B)

<caption>Fig. 20. Fiber orientation a<sub>11</sub> (A) and a<sub>22</sub> (B) components at B: comparison between the uncoupled/coupled simulations and experimental data</caption>
The relative error e between simulated and measured fiber orientations can be quantified by the following formula

$$e = \frac{\int_{-1}^{1} \|\mathbf{a}^{(\text{num})} - \mathbf{a}^{(\text{exp})}\| \mathrm{d}z}{\int_{-1}^{1} \|\mathbf{a}^{(\text{exp})}\| \mathrm{d}z}, \quad \|\mathbf{a}\| = \left( \sum_{(i,j) \in \mathcal{J} | a_{ij} |_{\mathrm{p}} \right)^{\frac{1}{p}}, \qquad (12)$$

where $p = 1$ and $\mathcal{J} = \{(1, 1), (2, 2)\}$ denotes the fiber orientation components used for error computation. Since both the experimental and numerical results are discrete in nature, they are linearly interpolated between two successive data points and then integrated in the z (gapwise) direction according to Eq. 12. The relative discrepancy as evaluated by Eq. 12 for the uncoupled and coupled simulation results is summarized in Table 2. The improvement introduced by a flow-fiber coupling viscosity is confirmed.

Finally, it is interesting to illustrate the consequences of such coupled simulations in overall mechanical performance predictions. Based on the micromechanical properties (resin and fiber) given in Tseng et al. (2017b) and the fiber orientation tensor calculated above, the overall elastic moduli of the composite can be predicted by the homogenization model presented in Tandon and Weng (1984) followed by the standard orientation averaging process (Advani and Tucker III, 1987). Since the directly measured elastic properties are not available, the experimental fiber orientation values $a_{11}$, $a_{22}$ and $a_{33}$ are used, by assuming that the off-diagonal components are zero. The orientation averaging is performed by the same ORT closure (Verweyst, 1998) used in the optimal scalar viscosity computation (Eq. 4). The homogenized Young's modulus in the 1 direction $E_{11}$ is presented in Fig. 21. Its average value in the

<table><caption>Table 2. Relative discrepancy e between numerical (uncoupled and coupled simulations) and experimental fiber orientation tensors at B</caption><thead><tr><th></th><th>Uncoupled simulation</th><th>Coupled simulation</th></tr></thead><tbody><tr><td>Relative discrepancy (%)</td><td>28</td><td>20</td></tr></tbody></table>

<figure><img src="image_31.png" /><figcaption>Fig. 21. Homogenized Young's modulus $E_{11}$ at B: comparison between the uncoupled/coupled simulations and experimental data</figcaption></figure>

thickness z direction is also summarized in Table 3. It can be observed that due to a slowed-down fiber alignment, the predicted tensile modulus is also decreased when coupling is considered. As expected, flow-fiber coupled simulations may also provide a correction to mechanical property predictions that could be overestimated in standard simulations due to an overaligned fiber orientation. Subsequent warpage simulation could hence benefit from this correction.

## 5 Conclusions and Future Work

From the point of view of rheological modeling of fiber-reinforced suspensions, this paper proposed a modified yet optimal version (Eq. 4) of the original viscosity model (Eq. 2) that reduces an anisotropic 4th order tensor to a scalar, even though it was initially motivated by computational implementation purposes. According to our preliminary theoretical analyses in section 2.2 and coupled simulation results presented in section 3, the optimal scalar model (Eq. 4) performs equally well compared to the original one and successfully captures elongational flows and fiber orientation components in the dominating shear plane.

Via these fiber-dependent viscosity models, flow kinematics and fiber evolutions are now tightly two-way coupled through rheological properties. In particular, a slowed down fiber alignment in the flow direction with a widened core region can be obtained, depending directly on the particular local deformation conditions and fiber states. The simulation results obtained in section 3 share qualitative similarities with several previous works referenced therein that also considered an orientation-dependent viscosity. Using the optimal scalar flow-fiber coupled model (Eq. 4), it is expected that the extensional flows will be better taken into account especially in the core region. Fiber orientation prediction for geometrically complex chunky parts (Kleindel et al., 2015) could be improved. In any case, the preliminary confrontation with experimental data presented in section 4 confirms the potential of such an orientation-dependent viscosity model.

In this paper, the flow-fiber coupling effects are evaluated using three different fiber orientation models available in Moldflow. Since once established the coupling works both ways, different fiber models may deliver different contributions to overall flow-fiber coupling. Attention must be taken when comparing results of coupled simulations with different orientation models. According to our results in section 3, it seems that the RSC and MRD models result in more coupling effects than the classical Folgar-Tucker model. Future work will be devoted to the interaction between fiber orientation models and fiber-dependent viscosity ones.

This paper reports some of the preliminary results of our research on fiber-induced viscosity modeling. It will be followed

<table><caption>Table 3. Through-thickness average $E_{11}$ (GPa) at B</caption><thead><tr><th>Uncoupled simulation</th><th>Coupled simulation</th><th>Experiment</th></tr></thead><tbody><tr><td>15.5</td><td>14.4</td><td>12.4</td></tr></tbody></table>
<figure><img src="image_32.png" /><figcaption>Fig. 22. Coupling between flow equations and fibers through fiber-induced viscosity models</figcaption></figure>

by additional communications that may include the following interesting theoretical and practical aspects:

* Comprehensive theoretical rheological analyses and computational details of the optimal scalar viscosity model (Eq. 4): further comparison with the original model in terms of approximation errors, implementation details in the Moldflow Insight API framework, computational efficiency ...

* Accounting for other additional fiber-induced rheological contributions. Different coupling may exist between flow equations and fiber evolution models, as indicated in Fig. 22. In this paper only orientation effects are considered. Future work could be for instance devoted to the coupling effects introduced by a non-homogeneous volume fraction distribution, which itself is predicted by a fiber migration model such as that proposed in Morris (2009).

* Implications of coupled injection molding simulations in subsequent integrative mechanical simulations that take into account the previous processing results: how the stress distribution is varied by using a fiber-flow coupling injection simulation.

* Further experimental validation of the proposed flow-fiber coupled simulations on structural parts. In section 4, the preliminary confrontation with experimental data yields promising results. However it was performed without considering packing effects and by using only the default RSC parameters. Future work will be devoted to these points in order to further decrease the discrepancy between numerical and experimental results. Optimization algorithms will be deployed to adjust and fine-tune the parameters of the fiber orientation models and the newly proposed coupling viscosity.

## References

Adam, L., Assaker, R., "Integrated Nonlinear Multi-Scale Material Modelling of Fiber Reinforced Plastics with Digimat: Application to Short and Continuous Fiber Composites", 11th World Congress on Computational Mechanics, Barcelona, Spain, p. 20–25 (2014)

Advani, S. G., Tucker III, C. L., "The Use of Tensors to Describe and Predict Fiber Orientation in Short Fiber Composites", J. Rheol., 31, 751–784 (1987), DOI:10.1122/1.549945

Akay, M., Barkley, D., "Flow-Aberrations and Weld Lines in Glass-Fibre Reinforced Thermoplastic Injection Mouldings", Plast. Rubber Compos. Process. Appl., 20, 137–149 (1993)

Arif, M. F., Saintier, N., Meraghni, F., Fitoussi, J., Chemisky, Y. and Robert, G., "Multiscale Fatigue Damage Characterization in Short Glass Fiber Reinforced Polyamide-66", Composites Part B, 61, 55–65 (2014), DOI:10.1016/j.compositesb.2014.01.019

Autodesk Moldflow Insight 2017 R2: Fiber Orientation Accuracy Validation Report, Autodesk Australia Pty Ltd., Kilsyth, VIC, Australia (2016)

Chung, S. T., Kwon, T. H., "Coupled Analysis of Injection Molding Filling and Fiber Orientation, Including In-Plane Velocity Gradient Effect", Polym. Compos., 17, 859–872 (1996), DOI:10.1002/pc.10679

Cintra Jr., J. S., Tucker III, C. L., "Orthotropic Closure Approximations for Flow-Induced Fiber Orientation", J. Rheol., 39, 1095–1122 (1995), DOI:10.1122/1.550630

Costa, F., Cook, P. and Pickett, D., "A Framework for Viscosity Model Research in Injection Molding Simulation, Including Pressure and Fiber Orientation Dependence", SPE ANTEC Conference, Orlando, FL, USA (2015)

Ericksen, J. L., "Anisotropic Fluids", Arch. Ration. Mech. Anal., 4, 231–237 (1959)

Folgar, F., Tucker III, C. L., "Orientation Behavior of Fibers in Concentrated Suspensions", J. Reinf. Plast. Compos., 3, 98–119 (1984), DOI:10.1177/073168448400300201

Jørgensen, J. K., Andreassen, E. and Salaberger, D., "The Effect of Fiber Concentration on Fiber Orientation in Injection Molded Film Gated Rectangular Plates", Polym. Compos., (2017), DOI:10.1002/pc.24698

Kleinindel, S., Salaberger, D., Eder, R., Schretter, H. and Hochenauer, C., "Prediction and Validation of Short Fiber Orientation in a Complex Injection Molded Part with Chunky Geometry", Int. Polym. Proc., 30, 366–380 (2015), DOI:10.3139/217.3047

Laun, H. M., "Orientation Effects and Rheology of Short Glass Fiber-Reinforced Thermoplastics", Colloid. Polym. Sci., 262, 257–269 (1984), DOI:10.1007/BF01410464

Laure, P., Silva, L. and Vincent, M., "Chapter 20 Modelling Short Fibre Polymer Reinforcements for Composites", in Composite Reinforcements for Optimum Performance, Boisse, P. (Ed.), Woodhead Publishing, Cambridge, UK (2011)

Lipscomb, G. G., Denn, M. M., Hur, D. U. and Boger, D. V., "The Flow of Fiber Suspensions in Complex Geometries", J. Non-Newtonian Fluid Mech., 26, 297–325 (1988), DOI:10.1016/0377-0257(88)80023-5

Mazahir, S., Vélez-García, G., Wapperom, P. and Baird, D., "Evolution of Fibre Orientation in Radial Direction in a Center-Gated Disk: Experiments and Simulation", Composites Part A, 51, 108–117 (2013), DOI:10.1016/j.compositesa.2013.04.008

Morris, J. F., "A Review of Microstructure in Concentrated Suspensions and Its Implications for Rheology and Bulk Flow", Rheol. Acta, 48, 909–923 (2009), DOI:10.1007/s00397-009-0352-1

Papathanasiou, T. D., "Chapter 4 Flow-Induced Alignment in Injection Molding of Fiber-Reinforced Polymer Composites", in Flow-induced Alignment in Composite Materials, Papathanasiou T. D., Guell, D. C. (Eds.), Woodhead Publishing, Cambridge, UK (1997)

Petrie, C. J. S., "Extensional Viscosity: A Critical Discussion", J. Non-Newtonian Fluid Mech., 137, 15–23 (2006), DOI:10.1016/j.jnnfm.2006.01.011

Phan-Thien, N., Graham, A., "A New Constitutive Model for Fibre Suspensions: Flow Past a Sphere", Rheol. Acta, 30, 44–57 (1991), DOI:10.1007/BF00366793

Phan-Thien, N., Zheng, R., "Chapter 3 Macroscopic Modelling of the Evolution of a Fibre Orientation during Flow", in Flow-Induced Alignment in Composite Materials, Papathanasiou, T. D., Guell, D. C. (Eds.), Woodhead Publishing, Cambridge, UK (1997)

Ranganathan, S., Advani, S. G., "A Simultaneous Solution for Flow and Fiber Orientation in Axisymmetric Diverging Radial Flow", J. Non-Newtonian Fluid Mech., 47, 107–136 (1993), DOI:10.1016/0377-0257(93)80047-F
Redjeb, A., Laure, P., Silva, L., Vincent, M. and Coupez, T., "Numerical Simulation of Fiber Orientation in Injection Molding Process", 21st Annual Meeting of the Polymer Processing Society, Leipzig, Germany (2005)

Rolland, H., Saintier, N. and Robert, G., "Damage Mechanisms in Short Glass Fibre Reinforced Thermoplastic during in situ Microtography Tensile Tests", Composites Part B, **90**, 365–377 (2016), DOI:10.1016/j.compositesb.2015.12.021

Tandon, G. P., Weng, G. J., "The Effect of Aspect Ratio of Inclusions on the Elastic Properties of Unidirectionally Aligned Composites", Polym. Compos., **5**, 327–333 (1984), DOI:10.1002/pc.750050413

Tseng, H.-C., Chang, R.-Y. and Hsu, C.-H., "Improved Fiber Orientation Predictions for Injection Molded Fiber Composites", Composites Part A, **99**, 65–75 (2017a), DOI:10.1016/j.compositesa.2017.04.004

Tseng, H.-C., Chang, R.-Y. and Hsu, C.-H., "Numerical Prediction of Fiber Orientation and Mechanical Performance for Short/Long Glass and Carbon Fiber-Reinforced Composites", Compos. Sci. Technol., **144**, 51–56 (2017b), DOI:10.1016/j.compscitech.2017.02.020

Tseng, H.-C., Chang, R.-Y. and Hsu, C.-H., "Numerical Predictions of Fiber Orientation for Injection Molded Rectangle Plate and Tensile Bar with Experimental Validations", Int. Polym. Proc., **33**, 96–105 (2018), DOI:10.3139/217.3404

Tucker III, C. L., "Flow Regimes for Fiber Suspensions in Narrow Gaps", J. Non-Newtonian Fluid Mech., **39**, 239–268 (1991), DOI:10.1016/0377-0257(91)80017-E

Verweyst, B. E., "Numerical Predictions of Flow-Induced Fiber Orientation in Three-Dimensional Geometries", PhD Thesis, University of Illinois, Urbana-Champaign (1998)

Verweyst, B. E., Tucker III, C. L., "Fiber Suspensions in Complex Geometries: Flow/Orientations Coupling", Can. J. Chem. Eng., **80**, 1093–1106 (2002), DOI:10.1002/cjce.5450800611

Vincent, M., Giroud, T., Clarke, A. and Eberhardt, C., "Description and Modeling of Fiber Orientation in Injection Molding of Fiber Reinforced Thermoplastics", Polymer, **46**, 6719–6725 (2005), DOI:10.1016/j.polymer.2005.05.026

Wang, J., O’Gara, J. F. and Tucker III, C. L., "An Objective Model for Slow Orientation Kinetics in Concentrated Fiber Suspensions: Theory and Rheological Evidence", J. Rheol., **52**, 1179–1200 (2008), PMid:18755969; DOI:10.1122/1.2946437

Wedgewood, A., Zhang, Z., Sulmoni, M. and Kang, S., "Addressing Practical Challenges in Developing Digimat Material Laws", Digimat Users' Meeting, Berlin, Germany (2017), PMid:29433682

Acknowledgements

This research was funded by the Région Île-de-France through the THERMOFIP research program (Pôle MOV’EO Convention 1709965). The authors would like to express again their sincere gratitude to Autodesk Moldflow developers and in particular Dr. Franco Costa for his guidance and valuable comments. Also, the Moldflow API simulations presented in this paper could not be performed without the kind and generous help of Alain Benchissou from APLICIT. Finally, we sincerely thank the reviewers for their constructive suggestions.

Date received: April 29, 2018

Date accepted: June 03, 2018

Bibliography

DOI 10.3139/217.3706

Intern. Polymer Processing

XXXIV (2019) 2; page 158–171

© Carl Hanser Verlag GmbH & Co. KG

ISSN 0930-777X