**2. Theory**

We briefly describe the different modeling techniques for fiber orientation estimation of short fiber reinforced thermoplastics (SFRT) in this section. Firstly, we give an overview of the fiber orientation tensor. Then, the common macroscopic models in the continuum scale are depicted starting from the basic models to the current state-of-the-art models. Finally, a mechanistic model based on the discrete element method (DEM) for the simulation on the particle level is discussed.

#### *2.1. Fiber Orientation Tensor*

A fiber can be represented as a unit vector under the assumption that the fibers are rigid, cylindrical, and have a uniform length and diameter [23]. Considering the above assumptions, we can define the fiber in space by a unit vector **p** and two angles (*θ*, *φ*) as shown in Figure 2.

**Figure 2.** Unit vector representing the orientation of a single rigid fiber.

The components of **p** can be represented in terms of *θ* and *φ* as:

$$p\_1 = \sin \theta \cos \phi \tag{1}$$

$$p\_2 = \sin\theta\sin\phi\tag{2}$$

$$p\_3 = \cos \theta.\tag{3}$$

A probability distribution function (PDF) *ψ*(**<sup>p</sup>**, *t*) can be used to represent the orientation of a population of fibers in space and time [9]. The PDF *ψ*(**<sup>p</sup>**, *<sup>t</sup>*)d**p** gives the probability of a fiber being directed between **p** and **p**+d**p** at a time *t*. The properties of the PDF are:

$$\mathcal{B}(\psi) = [0, 1] \tag{4}$$

$$
\oint \psi(\mathbf{p}, t) d\mathbf{p} = 1 \tag{5}
$$

$$
\psi(\mathbf{p}) = \psi(-\mathbf{p})\tag{6}
$$

$$\frac{D\psi}{Dt} = -\nabla\_S \cdot (\psi \dot{\mathbf{p}})\_\prime \tag{7}$$

where B indicates the image of the PDF. Equation (6) is only valid when the fibers are axis symmetrical and when neither of the ends of the fiber have a preferred head. Equation (7) is the continuity condition where ∇*S* is the gradient on the surface of the sphere formed by all possible orientations of an axis-symmetrical filler.

The representation of fiber orientation using the PDF is complex and inefficient for numerical simulation. A reduced representation was proposed for better computational efficiency [23]. This representation considers the moments of the PDF, which are calculated in the following way:

$$a\_{ij} = \int p\_i p\_j \psi(\mathbf{p}) \,\mathrm{d}\mathbf{p} \tag{8}$$

$$a\_{ijkl} = \int\_{\epsilon} p\_i p\_j p\_k p\_l \psi(\mathbf{p}) d\mathbf{p} \tag{9}$$

$$a\_{i\dots n} = \int p\_i \dots p\_n \psi(\mathbf{p}) \,\mathrm{d}\mathbf{p}.\tag{10}$$

They are called fiber orientation tensors. Since the PDF is of even order, as depicted in Equation (6), the odd ordered tensor integrals are equal to zero. Hence, we only consider the tensors of even order. There can be an infinite number of even ordered tensors. Among these even ordered tensors, mostly the second and fourth order tensors are used [23]. We introduce a compact notation **A** = *aij* and A = *aijkl* for further discussion.

The important properties of the orientation tensors are discussed only for the second and fourth order tensors in this work. The second order tensor is symmetric, and its trace is one [23],

$$\mathbf{A}\_{ij} = \mathbf{A}\_{ji} \tag{11}$$

$$\text{tr}\mathbf{A} = 1.\tag{12}$$

The fourth order tensor is symmetric for any pair of indices:

$$\mathbf{A}\_{ijkl} = \mathbf{A}\_{jikl} = \mathbf{A}\_{ijlk} = \mathbf{A}\_{kjil} = \mathbf{A}\_{ljki} = \mathbf{A}\_{ikjl} = \mathbf{A}\_{ilkj} \tag{13}$$

and it includes the entire information contained in the second order tensor [23]:

$$\mathbf{A}\_{ijkk} = \mathbf{A}\_{ij}.\tag{14}$$

The major advantages of using orientation tensors for computational purposes are that the unit sphere around a fiber does not need discretization and the three-dimensional calculations can be feasibly done independently of the basis system [23].

#### *2.2. Continuum-Based Models for Fiber Orientation*

A substantial amount of research work has been done to model the fiber orientation for fiber reinforced polymers since they play an important role in determining the structural strength of final parts. One of the very first models depicting the motion of a single fiber in a non-turbulent Newtonian fluid was given by Jeffrey [2]. The fiber is modeled as a rigid ellipsoid, which does not exhibit bending, nor attrition.

$$
\dot{\mathbf{p}} = \mathbf{W} \cdot \mathbf{p} + \zeta (\mathbf{D} \cdot \mathbf{p} - \mathbf{D} : \mathbf{pp} \mathbf{p}) \tag{15}
$$

$$\mathbf{W} = \frac{1}{2} (\nabla \mathbf{u} - \nabla \mathbf{u}^T) \tag{16}$$

$$\mathbf{D} = \frac{1}{2} (\nabla \mathbf{u} + \nabla \mathbf{u}^T)\_\prime \tag{17}$$

where **W** is the vorticity tensor and **D** is the rate of strain tensor with ∇**u** being the flow gradient. The term *ξ* represents the particle shape parameter where *ξ* = *λ*2−1 1+*λ*<sup>2</sup> and *λ* is the fiber aspect ratio.

The Jeffrey model is valid only in dilute solutions, where the dominant mechanism is fluid-fiber interaction since it was developed as a model for single fiber motion [2]. Industrial thermoplastics are typically highly charged in fiber (concentrated regime), where fiber-fiber interaction is dominant. Hence, newer models taking into account the fiber-fiber interactions were developed. All models for the prediction of fiber orientation, except the Jeffrey model, are phenomenological.

The Folgar and Tucker (FT) model [9] uses an additional diffusion term to take into account the fiber-fiber interaction in semi-dilute and concentrated regimes. This model considers the fibers as rigid and uniform cylinders, but large enough to prevent Brownian motion. The fluid is incompressible and viscous enough in order to neglect the inertia of particle and its buoyancy. This model was developed in terms of the rate of change of the second order oriented tensor by Advani and Tucker [23]:

$$\mathbf{^0\frac{D}{Dt}} = \mathbf{\dot{A}}\; \mathbf{\dot{A}}\; = \mathbf{\dot{A}}^h + \mathbf{\dot{A}}^d\tag{18}$$

$$\mathbf{A}^{\hbar} = \left(\mathbf{W} \cdot \mathbf{A} - \mathbf{A} \cdot \mathbf{W}\right) + \zeta \left(\mathbf{D} \cdot \mathbf{A} + \mathbf{A} \cdot \mathbf{D} - 2\mathbf{A} : \mathbf{D}\right) \tag{19}$$

$$\mathbf{A}^{\mathbf{d}} = 2C\_{\mathbf{l}}\gamma(\mathbf{I} - 3\mathbf{A}),\tag{20}$$

where **A**˙ h corresponds to the hydrodynamic effects and **A**˙ d is the diffusion term, which takes care of the fiber-fiber interactions. The term *C*I is the interaction coefficient, and *γ*˙ = √2**D** : **D** is the magnitude of the rate of the strain tensor. We see that in Equation (19), the fourth order orientation tensor is used. Hence, a closure approximation is necessary. There are several types of closure approximations proposed in the literature as for example fitted closures such as orthotropic closures [24–26], exact closures [27–29], or simple closures [23,30]. In our work, we use a fitted closure, the invariant-based optimal fitted (IBOF) closure approximation [25] based on the findings of Kugler et al. [15].

The Folgar–Tucker model predicts a faster evolution of fiber orientation when compared to experimental results [3]. Therefore, newer models introducing a retardation of the fiber orientation evolution have been developed. All the aforementioned assumptions are valid for the following models as they are based on [9].

The reduced strain closure (RSC) model introduced by Wang et al. [3] is based on the spectral decomposition of **A** where **A** = ∑3 *i*=1 *λi***e***i***e***i*. The modified growth rate is given as:

$$
\lambda\_i^{\mathsf{RSC}} = \kappa \lambda\_i \tag{21}
$$

$$
\dot{\mathbf{e}}\_i^{\text{RSC}} = \dot{\mathbf{e}}\_{i\prime} \tag{22}
$$

where *κ* is a constant such that *κ* ∈ [0, 1]. The Folgar–Tucker Equation (18) along with the modified growth rate Equation (21) can be expressed as:

$$\mathbf{A} = \mathbf{A}^{\text{RSC}} + \kappa \mathbf{A}^d \tag{23}$$

$$\mathbf{A}^{\mathsf{RSC}} = \mathbf{W} \cdot \mathbf{A} - \mathbf{A} \cdot \mathbf{W} + \mathfrak{J} \left[ \mathbf{D} \cdot \mathbf{A} + \mathbf{A} \cdot \mathbf{D} - 2 (\mathbf{A} + (1 - \kappa)(\mathbf{L} - \mathbf{M} : \mathbf{A})) : \mathbf{D} \right] \tag{24}$$

$$\mathbf{L} = \sum\_{i=1}^{3} \lambda\_i \mathbf{e}\_i \mathbf{e}\_i \mathbf{e}\_i \mathbf{e}\_i \tag{25}$$

$$M = \sum\_{i=1}^{3} \mathbf{e}\_i \mathbf{e}\_i \mathbf{e}\_i \mathbf{e}\_i. \tag{26}$$

The output of the Folgar–Tucker model could not be properly fitted with the experimental data [15]. Hence, Phelps and Tucker [6] formulated an anisotropic rotary diffusion model applicable to both short and long fiber thermoplastics. In this model, the interaction coefficient *CI* is substituted by a rotary diffusion tensor **C**. It can be additionally noted that in this model, the state of fiber orientation has an influence on the rotary diffusion effect.

$$\mathbf{A} = \mathbf{A}^h + \mathbf{A}^{\text{ARD}} \tag{27}$$

$$\mathbf{A}^{\rm ARD} = \dot{\gamma} [2\mathbf{C} - 2\text{tr}(\mathbf{C})\mathbf{A} - 5(\mathbf{C} \cdot \mathbf{A} + \mathbf{A} \cdot \mathbf{C}) + 10\mathbf{A} : \mathbf{C}].\tag{28}$$

To reduce the amount of phenomenological parameters, the pARD model was developed by Tseng et al. [7]. In their model, the rotary diffusion tensor **C** is represented as:

$$\mathbf{C} = \mathbf{C}\_{\text{I}} \mathbf{R}\_{\text{A}} \begin{pmatrix} D\_1 & 0 & 0 \\ 0 & D\_2 & 0 \\ 0 & 0 & D\_3 \end{pmatrix} \mathbf{R}\_{\text{A}'}^\top \tag{29}$$

where **R**A is the eigen matrix. The work of [7] used *D*1 = 1, *D*2 = c and *D*3 = 1 − c to reduce the number of parameters. This implies that *C*I determines the rotary diffusion factor in the first principal direction for fiber orientation. The second and the third principal fiber orientation directions are scaled with the factors c · *C*I and (1 − c) · *C*I, respectively.

An ARD-RSC model is also introduced based on the work of [6] to further retard the kinetics of the ARD model:

$$\begin{split} \mathbf{\dot{A}}^{\text{ARD},\text{RSC}} = & \mathbf{W} \cdot \mathbf{A} - \mathbf{A} \cdot \mathbf{W} + \zeta [\mathbf{D} \cdot \mathbf{A} + \mathbf{A} \cdot \mathbf{D} - 2(\mathbf{A} + (1 - \kappa)(\mathbb{L} - \mathbf{M} : \mathbb{A})) : \mathbf{D}] \\ &+ \dot{\gamma} [2(\mathbf{C} - (1 - \kappa)\mathbf{M} : \mathbf{C}) - 2\kappa(\text{tr}(\mathbf{C})) \mathbf{A} - 5(\mathbf{C} \cdot \mathbf{A} + \mathbf{A} \cdot \mathbf{C}) \\ &+ 10[\mathbf{A} + (1 - \kappa)(\mathbb{L} - \mathbf{M} : \mathbb{A})] : \mathbf{C}]. \end{split} \tag{30}$$

which reduces to the ARD model if *κ* = 1. Another novel approach is the MRD (Moldflow rotational diffusion) model from the work of Bakharev et al. [8],

$$\mathbf{A} = \mathbf{A}^h + \mathbf{A}^{\text{MRD}} \tag{31}$$

$$\mathbf{A}^{\text{MRDD}} = 2\gamma(\mathbf{C} - \text{tr}(\mathbf{C})\mathbf{A}).\tag{32}$$

#### *2.3. Particle-Based Mechanistic Model*

We use a mechanistic model for the determination of fiber positions under the influence of different flow types and varying boundary conditions. The fiber model is based on the discrete element method (DEM) and is derived from the works of Perez [1] and Lindström [31].

A fiber is represented by a chain of cylindrical rods. The cylindrical rods are rigid and interconnected via ball and socket joints as shown in Figure 3. The cylindrical rods are known as segments, and the ends of each segmen<sup>t</sup> are called nodes. In each segment, the positions **x***i*, velocities **u***i*, and angular velocities *ωi* are stored at every node *i*. The fibers are suspended in a constant viscosity fluid, which is submitted to an external flow field **U**<sup>∞</sup>. The system is partially coupled such that the hydrodynamic forces on fibers due to the flow are taken into account, but the force exerted on the fluid by the fibers is neglected. The different forces acting on each segmen<sup>t</sup> *i* are: hydrodynamic forces **F***Hi* from the fluid, the interaction force with neighboring segmen<sup>t</sup> *j* **F***Cij*, and intra-fiber force from the

adjacent segmen<sup>t</sup> *i* **X***i*. **M***bi* and **<sup>M</sup>***bi*+<sup>1</sup> are the bending moments. These forces and moments are used to formulate the translational equation of motion (33) and the rotational equation of motion (34).

$$\mathbf{F}\_i^H + \sum\_j \mathbf{F}\_{ij}^C + \mathbf{X}\_i - \mathbf{X}\_{i+1} = \mathbf{0} \tag{33}$$

$$\mathbf{T}\_i^H - \mathbf{r}\_i \times \mathbf{X}\_{i+1} + \sum\_j d\_{ij} \times \mathbf{F}\_{ij}^C + \mathbf{M}\_i^b - \mathbf{M}\_{i+1}^b = 0 \tag{34}$$

$$
\mathbf{u}\_i + \boldsymbol{\omega}\_i \times \mathbf{r}\_i - \mathbf{u}\_{i+1} = \mathbf{0},\tag{35}
$$

where **r***i* is the segmen<sup>t</sup> vector and *dij* is the minimum distance between two neighboring segments, as shown in Figure 3. Equation (35) is necessary to enforce connectivity if the fibers are divided into multiple segments.

**Figure 3.** Fiber represented as a chain of segments. A neighboring fiber is also represented (adapted from [31]).

The definitions used in this model for the hydrodynamic and bending forces are similar to [1]. The interaction force between fibers is modeled based on the work of [31]. The interaction between fibers is divided into three regimes: Coulomb friction and normal contact forces **F***mc ij* when a mechanical contact exists between fiber segments, transitional forces **F***trans ij* at distances lower than the fiber surface roughness *δsur*, and lubrication forces **F***lub ij* , which takes place at distances larger than *δsur*, but shorter than the threshold *δlub* (Equation (36)).

$$\mathbf{F}\_{ij}^{\mathbb{C}} = \begin{cases} \mathbf{F}\_{ij}^{\text{lab}} & \delta\_{\text{sur}} < d\_{ij} < \delta\_{\text{lub}} \\ \mathbf{F}\_{ij}^{\text{trans}} & \mid d\_{ij} \mid < \delta\_{\text{sur}} \\ \mathbf{F}\_{ij}^{\text{mc}} & d\_{ij} < -\delta\_{\text{sur}} \end{cases} \tag{36}$$

In the case of the inclusion of walls in the system, the interaction among fibers and walls is computed in the same manner as the computation of the interaction between fibers. The solution of Equations (33)–(35) gives the position of each fiber at every time step. In post-processing, from fiber positions at every time step, the second order orientation tensor can be computed as follows [32]:

$$\mathbf{A}\_{i\bar{j}}(t) = \frac{\sum\_{n=1}^{N} p\_{n,\bar{j}}(t) p\_{n,\bar{i}}(t)}{N},\tag{37}$$

where *N* is the total number of fibers in the system, *t* is time, and *pn*(*t*) is the position of the fiber *n* at time *t*. The evolution of the fiber orientation in the directions of the global system can be obtained from the diagonal components of the orientation tensor, **A**11, **A**22, and **A**33, which correspond to the fraction of total fibers oriented preferentially in the *X*-, *Y*-, and *Z*-directions, respectively.

#### **3. Fiber Orientation in Elongational Flows**

The mechanistic model described in the previous section is now used to simulate the effect of elongational flow on fiber orientation.

#### *3.1. Simulation Method and Comparison with the Experiment*

Kugler et al. [32] used the mechanistic model to simulate the effect of shear flow using periodic boundary conditions. In our work, this mechanistic model is extended for simulating also fiber orientation under elongational flow. The simulation is performed with the help of a representative volume element (RVE) created using the method developed by Schneider [33]. The elongational flow on the RVE is applied by emulating walls in two directions and imposing the following velocity field:

$$
\hbar u\_x = \dot{\epsilon} x, \; u\_y = -\dot{\epsilon} y, \; u\_z = 0. \tag{38}
$$

where ˙ is the rate of elongation. This flow field can be interpreted as a fluid flowing from the *y*-direction towards the *x*-direction, but constrained in the *z*-direction. Walls are defined on the top and bottom of the RVE, which are moving with velocities *vupper w* and *vlower w* towards the negative *y*-direction and the positive *y*-direction, respectively. They are given according to:

$$
v\_{w|y}^{upper} = -\frac{1}{2}\dot{\varepsilon}\,h(t),\ v\_{w|y}^{lower} = \frac{1}{2}\dot{\varepsilon}\,h(t),\tag{39}$$

where *h*(*t*) is the distance between the walls as they approach each other during the simulation. It should be noted that the wall velocities *vw*|*y* reduce as the height of the cell decreases. The wall velocities after a certain time become too small for the numerical procedure to continue. Hence, a lower limit is added to the wall velocities beyond which they are kept constant. This lower limit depends on the elongational rate and time. We expect this to have a negligible influence on the estimated fiber orientation, since only a small amount of fibers close to the walls are directly affected. The walls are also positioned 2 × 10−<sup>4</sup> m away from the edge fibers in the initial RVE configuration. In that way, we avoid any initial instability arising out of a fiber with ends outside the RVE, where the velocities cannot be numerically resolved.

Additionally, walls in the *xy*-plane are included. Although there is no flow defined in the *z*-direction, it prevents fibers from leaving the RVE. The RVE along with the boundary walls is shown in Figure 4.

**Figure 4.** Visualization of the boundary walls of the RVE for elongational flow. (**a**) Before the start of simulation. (**b**) After the end of simulation.

There have been no experimental works so far conducted for elongational flow using short fiber reinforced thermoplastics. However, there were two experiments conducted by Lambert et al. [17] and Lambert and Baird [18] using long fibers, under the influence of elongational flow. In the following sections, we use the mechanistic model for the simulation of elongational flow with long fibers to compare the results with the experimental work.

3.1.1. Validation with Experimental Non-Lubricated Squeeze Flow

Lambert et al. [17] used a commercial composite Verton MV006S (30% wt. glass fiber reinforced polypropylene (PP)) for the experiment. Samples were fabricated by compression molding; for a detailed description, refer to [17]. The squeeze flow experiment was done using a custom-built device described in [17]. There was no lubricant used at the interaction of the mold walls and the samples. The process settings for the experiment are given in Table 1.

**Table 1.** Process settings of the non-lubricated squeeze flow experiment for a 30% wt. glass fiber reinforced polypropylene (PP) in [17].


The mechanistic model defined in the previous section is now used to simulate the elongational flow with long fibers mimicking the settings in Table 2. The initial fiber orientation obtained from [17] was **A**11 = **A**33 = 0.45 and **A**22 ∼= 0.1 (averaged over three samples). This initial fiber orientation was used to generate an RVE with a squared cross-section of 4 mm<sup>2</sup> using the method given in [33]. The RVE cross-section was smaller than the experimental sample due to computational constraints. The viscosity was approximated as three times the zero shear viscosity of the pure PP matrix to account for elongational viscosity.

**Table 2.** Simulation parameters for comparison with the non-lubricated squeeze flow experiment.


The flow field and boundary conditions were similar to the ones we described above. It is important to note that during the simulation, we considered no friction between the walls of the RVE and the fibers, as opposed to the experimental setup. The RVE was divided into five layers over the thickness direction (*y*-direction) for post-processing. We plot in Figure 5 the diagonal components of the orientation tensor **A**11, **A**22, and **A**33 over the normalized thickness position obtained from the layers with the experimental results found in [17].

The simulation results showed good agreemen<sup>t</sup> in the middle of the thickness. However, they slightly under predicted **A**33 at the core. The overall shape of the fiber orientation profile from the experiment was not well captured by the simulation. This was evident from the mismatch of the simulation and experimental results at the boundaries. This can be attributed to the fact that we did not account for friction at walls in the simulation.

**Figure 5.** Final fiber orientation through thickness after a elongation time of 2 s at a constant strain rate of 0.5 s<sup>−</sup><sup>1</sup> with experimental data from [17] in points and simulation data in continuous lines.

#### 3.1.2. Validation with Experimental Constant Planar Extension

Lambert and Baird [18] used a 10% wt. glass fiber reinforced PP for the constant planar extension experiment. All samples were made using compression molding. During the course of the experiment, a silicone lubricant (DCF 203, Dow Corning, Midland, MI, USA) with low viscosity was applied to prevent friction and reduce the shear between the mold wall and the melt [18]. The samples were squeezed in a custom-built device at a constant strain rate. The process settings are described in Table 3, and the experiment is explained in detail in [18].

**Table 3.** Process settings of the constant planar extensional flow experiment in [18].


The introduced mechanistic model was then used for comparison with the experimental results from [18]. Due to computational constraints, the RVE used had a reduced squared cross-section of 9 mm2. The simulation parameters are listed in Table 4. The RVE was compressed in the *y*-direction at a constant strain rate during 40 s.

**Table 4.** Simulation parameters for comparison with the planar extensional flow experiment.


The second order fiber orientation tensor **A** was computed at each time step. The evolution of the fiber orientation in the principal directions **A**11, **A**22, and **A**33 is plotted against time along with the experimental results from [18] in Figure 6.

**Figure 6.** Fiber orientation evolution under constant planar extension at a strain rate of 0.05 s<sup>−</sup><sup>1</sup> with experimental results from [18] in points and simulation results in continuous lines.

The good agreemen<sup>t</sup> of the simulation with the experimental results for constant planar extensional flow can be seen. The evolution of fiber orientation occurs at a very fast rate for this type of flow, which was confirmed by both the experimental and the mechanistic model results in our case. Since there was no parameter fitting involved for the output from our mechanistic model, it can be said that it gave a good estimation of fiber orientation under elongational flow. This model can now be used to predict fiber orientation evolution in elongational flows for various different materials.

#### *3.2. Numerical Analysis: Factors Influencing Fiber Orientation in Elongational Flow*

In this section, we study the effects of flow rate, several fiber descriptors, and polymer viscosity on the fiber orientation evolution in elongational flow.

To study the effect of the elongational rate, the simulation was performed with the following elongational rates: 1.0 s<sup>−</sup><sup>1</sup> and 100.0 s<sup>−</sup>1. Plotted over strain, there was almost no difference in the final fiber orientation for the analyzed flow rates, as seen from Figure 7. Thus, it seemed that the elongational rate scaled linearly with orientation evolution.

**Figure 7.** Fiber orientation evolution under elongational flow with different elongation rates.

We applied a similar procedure for the evaluation of the influence of fiber length. The fiber radius, the fiber volume fraction, and fluid viscosity were kept constant, and only the fiber length was varied. The different fiber lengths ranging from short fibers to long fibers were 2.5 × 10−<sup>4</sup> m, 8 × 10−<sup>4</sup> m, and 1.5 × 10−<sup>3</sup> m. It can be observed from Figure 8 that the fiber length had a negligible effect on the fiber orientation.

**Figure 8.** Fiber orientation evolution under elongational flow with different fiber lengths.

The influence of the volume fraction on fiber orientation was investigated using a fixed fiber length, fiber radius, and fluid viscosity. We chose short fibers having a length of 2.5 × 10−<sup>4</sup> m and a radius of 5 × 10−<sup>6</sup> m, since they are used extensively in industrial parts. The fluid viscosity was fixed at 900 Pa s. The RVE had a squared cross-section of 2.5 mm<sup>2</sup> and a height of 3 × 10−<sup>3</sup> m. The volume fraction of each RVE is listed in Table 5. An elongational rate of 0.05 s<sup>−</sup><sup>1</sup> was applied until reaching a strain value of two. The fiber orientation evolution over strain is depicted in Figure 9.

**Table 5.** Volume fractions and number of fibers in the respective RVEs for analyzing the influence of fiber content on fiber orientation under elongational flow.

**Figure 9.** Fiber orientation evolution under elongational flow with different RVEs having different volume fractions.

There was no remarkable difference in fiber orientation evolution with different volume fractions, as seen from Figure 9. In the case of close examination of the evolution curves, we saw that at higher volume fractions of 30% and 25%, the fiber orientation evolution was slightly faster in the principal directions **A**11 and **A**33.

As the last factor, the effect of the polymer matrix viscosity on the fiber orientation evolution was studied. A flow rate of 0.5 s<sup>−</sup>1, a volume fraction of 18%, a fiber radius of 5 × 10−<sup>6</sup> m, and a fiber length of 2.5 × 10−<sup>4</sup> m were used. The RVE height was 3 × 10−<sup>3</sup> m. The viscosity was varied from 300 Pa s to 1200 Pa s. The fiber orientation evolution with respect to strain is plotted in Figure 10.

**Figure 10.** Fiber orientation evolution under elongational flow with different matrix viscosities.

We can observe from Figure 10 that there was a difference in fiber orientation evolution between the lower viscosities (300 Pa s and 600 Pa s) and the higher viscosities (900 Pa s and 1200 Pa s). A possible explanation is that the entire system became less stiff as the fluid viscosity reduced. This was evident given that a lower viscosity reduced the lubrication forces arising in the system, and hence, the fibers could orient more easily under the influence of flow.

In summary, we observed that there was a negligible effect of elongational rate and fiber length. Fiber volume fraction and matrix viscosity had a slight effect on the fiber orientation evolution under elongational flow, but this variation was minor when compared to the fiber orientation evolution under other types of flow (e.g., simple shear). Therefore, it seems that the fiber orientation evolution under elongational flow is universal and independent of strain rate intensity and material descriptors.

It has to be clarified that all effects were studied separately. Additional research is necessary to determine the effect of volume fraction variations with long fibers as well as the effect of fiber length at high volume fractions. The combination of high volume fraction and high fiber length has not been studied, since fibers have been modeled as rigid. The stated combination seems critical under this assumption. Overall, the assumption of rigid fibers seems less severe in elongational flow than in shear flow, since fibers do not perform orbits. In shear flow, the critical length where fiber flexibility has to be incorporated for single fibers movement is around an aspect ratio of 100 [34]. In elongational flow, we assume that the movement of single fibers can be modeled as rigid. Nevertheless, we do not assume that this assumption holds true with increased fiber contact at high volume fractions for long fibers.

#### *3.3. Evaluation of Macroscopic Fiber Orientation Models in Elongational Flows*

We conducted a calibration of macroscopic fiber orientation models in elongational flows. Since the mechanistic model is able to reproduce the fiber orientation phenomenon under elongational flow, a virtual elongational test for a 30 wt.% glass fiber reinforced PP was conducted. The simulation settings are stated in Table 6.

**Table 6.** Simulation settings for the virtual elongational test for PP-GF30.


A generic algorithm was used to determine the best fit of the diverse macroscopic fiber orientation model parameters. It can be observed that a scalar diffusion rate was sufficient to fit the evolution data in elongational flows. Figure 11 displays a comparison of the best fits of the FT, RSC, and pARD-RSC model. The other models are omitted since they show similar results. Consequently, only one parameter was fitted *CI* = 0.0019. In all models, the parameters were set such that the FT model was retrieved, e.g., *κ* = 1 or *D*2 = *D*3 = 1.

**Figure 11.** Comparison of the best fit of the FT, RSC, and pARD-RSC model with the virtual elongational test for PP-GF30. Simulation settings are stated in Table 6.

The fit is repeated for the different volume fractions and viscosity data in Section 3.2 to evaluate whether these changes at the micro scales had an impact on the macroscopic fiber orientation model parameters.

For low viscosities (300 Pa s to 600 Pa s), the fitted FT model slightly underpredicted the evolution speed of the orientation. The interaction coefficient for all viscosities varied between 0.001 and 0.003. Thus, we assumed that one constant fiber orientation model parameter set under elongational flow was valid for a broad variety of materials.

#### **4. Influence of Flow Type on Fiber Orientation**

In this section, we propose a macroscopic model for incorporating the flow-type effect on fiber orientation and implement it in the injection molding software Autodesk Moldflow Insight Scandium® 2019. The existing and newly implemented models from the work of Kugler et al. [15] took into account the effect of shear flow on fiber orientation. In this work, we introduce a novel model approach to consider the effect of both shear and elongational flows for the estimation of fiber orientation in an injection-molded part.

#### *4.1. Flow-Dependent Fiber Orientation Model Definition*

During injection molding, shear and elongational flows occur. To determine the main flow regime of a velocity field at a given position and time, an objective measure is necessary. We used the Manas-Zloczower number Mz [35]. Instead of the Manas-Zloczower number, the objective measure *β* by Astarita [36] can also be used.

The Manas-Zloczower number is defined by:

$$\mathbf{M}\_{\mathbf{z}} = \frac{||\mathbf{D}||}{||\mathbf{D}|| + ||\mathbf{W}||} \tag{40}$$

where **D** is the magnitude of the rate of the strain tensor and **W** is the magnitude of the vorticity tensor. The magnitudes are calculated by using the square roots of the second invariants of the rate of strain and vorticity tensor. Consequently, the Manas-Zloczower number Mz is objective. The second invariant of the rate of strain tensor is defined by tr **D**<sup>2</sup> [37], so:

$$\|\|\mathbf{D}\|\| = \sqrt{\text{tr}\,\mathbf{D}^2}.\tag{41}$$

The second invariant of the vorticity tensor is defined as tr **W**<sup>2</sup> = 12 **Ω** 2, where **Ω** = ∇**u** is the vorticity vector [37]. Therefore, the magnitude of the vorticity tensor is defined by:

$$\|\mathbf{W}\| = \sqrt{\frac{1}{2} \|\mathbf{\Omega}\|^2}. \tag{42}$$

The Manas-Zloczower number characterizes the flow type in the following way:

$$\mathsf{M}\_{\mathsf{Z}} = \begin{cases} 1 & \text{elongational flow} \\ 0.5 & \text{shear flow} \\ 0 & \text{rotational flow} \end{cases} \tag{43}$$

In our work, we consider the effect of shear and elongational flow. Hence, we applied a lower limit of 0.5 to the flow type descriptor M<sup>z</sup>. If at a given point during the flow, Mz < 0.5 holds true, then it is set to 0.5, assuming that the flow is controlled by shear. We believe that this hypothesis is weak given the low frequency of rotational flow in front of shear and elongational flows in injection molding.

We assumed that fiber orientation for mixed shear and elongational flow is a linear combination of the shear- and elongational-driven fiber orientation models. This leads to the following flow-dependent fiber orientation model:

$$\mathbf{A} = 2(\mathbf{M}\_{\mathbf{z}} - 0.5)\mathbf{A}\_{\mathbf{c}} + 2(1 - \mathbf{M}\_{\mathbf{z}})\mathbf{A}\_{\mathbf{s}}.\tag{44}$$

In pure shear, this equals **A** ˙ = **A** ˙ *s*, and in pure elongational flow, **A** ˙ = **A** ˙ *e*. The derivatives **A** ˙ *s* and **A** ˙ *e* can be calculated using a macroscopic fiber orientation model with the respective parameters for shear and elongation. Based on the findings of Kugler et al. [15], we used the pARD-RSC model to model fiber orientation in shear flow. From the mechanistic model simulation and the respective macroscopic fiber orientation model fits in Section 3.3, it was shown that the FT model is sufficient in elongational flows. For the sake of consistency, we also used the pARD-RSC model in elongation. The number of parameters for the pARD-RSC model is then reduced in the case of elongational flow by setting *D*1 = *D*2 = *D*3 = 1 in Equation (29). Additionally, no retarding rate is necessary; therefore, we set *κ* = 1. Consequently, only the interaction coefficient *CI* has to be fitted.

This leads to the following definition of the flow-dependent fiber orientation model:

$$\dot{\mathbf{A}}^{\text{pARD}-\text{RSC}} = 2(\mathcal{M}\_{\text{Z}} - 0.5)\dot{\mathbf{A}}\_{\text{c}}^{\text{pARD}-\text{RSC}} + 2(1 - \mathcal{M}\_{\text{Z}})\dot{\mathbf{A}}\_{\text{s}}^{\text{pARD}-\text{RSC}}.\tag{45}$$

with:

$$\mathbf{A}\_{\varepsilon}^{\mathsf{p}\mathsf{A}\mathsf{R}\mathsf{D}-\mathsf{R}\mathsf{S}\mathsf{C}} = \mathbf{A}^{\mathsf{p}\mathsf{A}\mathsf{R}\mathsf{D}-\mathsf{R}\mathsf{S}\mathsf{C}}(\mathsf{C}\_{l}^{\mathsf{c}},1,1,1) \tag{46}$$

$$\mathbf{A}\_s^{\mathsf{P}\mathsf{A}\mathsf{R}\mathsf{D}-\mathsf{R}\mathsf{S}\mathsf{C}} = \mathbf{A}^{\mathsf{p}\mathsf{A}\mathsf{R}\mathsf{D}-\mathsf{R}\mathsf{S}\mathsf{C}}(\mathbf{C}\_1^s, \mathbf{D}\_2^s, \mathbf{D}\_3^s, \mathbf{x}^s). \tag{47}$$

This linearly interpolated fiber orientation model was implemented in Autodesk Moldflow Insight Scandium® 2019 using the Solver API feature. The flow-dependent model can estimate the fiber orientation of any fiber reinforced plastic part using user-defined model parameters.

#### *4.2. Workflow for Fiber Orientation Estimation Using the Flow-Dependent Model*

In this section, we propose a workflow for the identification of the shear and elongational parameters for the flow-dependent fiber orientation model.

The parameters for shear flow can either be experimentally fitted like in [15] or determined by a virtual shearing test with a mechanistic fiber simulation [38].

The parameters for elongational flow can be determined by the virtual elongational test proposed in Section 3. Additionally, the experimental setup of [17] can be used. This workflow is displayed in Figure 12.

**Figure 12.** Workflow for parameter identification of the flow-dependent fiber orientation model.

#### *4.3. Flow-Dependent Model Parameters for a PBT-GF30*

In this section, the workflow is applied to define the parameters for a PBT-GF30. In shear, we used the experimental parameters determined in [15]. In elongation, we used the parameters from the virtual elongational test with the mechanistic model from Section 3.3. The test was conducted for a PP-GF30, but since the variation in the polymer matrix viscosity was small, we supposed that the same parameter set under elongational flow was also valid for our PBT-GF30. The parameters of the flow-dependent fiber orientation model PBT-GF30 are given in Table 7.


#### **5. Results and Discussion**

In this section, we evaluate the accuracy of diverse fiber orientation models by comparing with experimental fiber orientation determined by microcomputed tomography (μCT) in different parts. We use the following macroscopic fiber orientation models:


We simulated a total of three geometries with the help of these models. All parts were injection-molded with a PBT-GF30, the same material referred to in Section 4.3. The parts included a 2 mm plate, a window gear lift housing, and a sensor cover. μCT scans were conducted at different positions in each geometry. The experimental fiber orientations were calculated from the analysis of the μCT data. The geometries and the positions of the μCT scans are depicted in Figure 13. The sensor cover cannot be fully displayed due to company restrictions. Position 3 at the sensor cover was located close to the gate. The locations were chosen in order to cover different flow lengths, different wall thicknesses, flow singularities like weld lines, and areas with specific mechanical requirements under assembling and/or operation.

The pARD-RSC model parameters for shear and elongational flow are given in Table 7. All geometries were simulated using a fill-and-pack analysis in Autodesk Moldflow Insight® 2019. The RSC and MRD models are pre-defined in Autodesk Moldflow®. The different versions of the pARD-RSC model were implemented in an extended technology preview of the software (Project Scandium) using the Solver API feature. For a validation of the implementation, refer to [15].

The simulation of the different geometries with the flow-dependent model requires the calculation of the Manas-Zloczower number during flow. We calculated the Manas-Zloczower number for each node at every time step during the filling of the mold. The coherence of the computation of the Manas-Zloczower number was evaluated with the help of the plate.

The Manas-Zloczower number was written as an output during the analysis using the Solver API. Figure 14a shows the cross-section of the plate with the different values of Mz after 0.41 s of filling. For information, the cavity was completely filled after 0.63 s. We could observe that Mz exhibited higher values at the core of the sprue and plate. The value of Mz was close to 0.5 near the boundaries of the plate. This agreed with the fact that shear flow was dominant at the boundaries due to the friction with the mold walls. The flow gradually changed to elongational flow towards the core of the plate. This can also be observed in Figure 14b where the Manas-Zloczower number is plotted against the thickness of the plate at Position 1 (Figure 13) once the node is fully filled for the first time.

**Figure 13.** Different geometries with the position of the μCT scans. (**a**) Plate. (**b**) Window lift drive housing. (**c**) Sensor cover (partly displayed).

**Figure 14.** Evaluation of the Manas-Zloczower number calculation. (**a**) Cross-section of plate after 0.41 s of filling time. (**b**) Graph of Mz against the thickness of the plate at Position 1 once the node is fully filled for the first time.

We evaluated the performance of the different models with the help of the root mean squared deviation (RMSD) between the eigenvalues of the fiber orientation tensor obtained from the μCT scans and the eigenvalues of the final fiber orientation tensor from simulation. The eigenvalues of **A** were used in the comparison instead of the diagonal entries **A***ii* for *i* = 1, 2, 3 due to the non-coinciding coordinate axes of the μCT analyses and those of the simulation. We depict the RMSD values of each model at all μCT scan positions and the average over all the positions in Table 8. The columns Shear and Elongation represent the pARD-RSC model with parameters from shear flow and elongational flow, respectively.

We can observe from the RMSD values for the different parts that the flow-dependent model gives similar or better results in comparison to the standard fiber orientation models in Moldflow. In the plate, models without a retarding rate yield better results (MRD and elongational pARD-RSC model). On the contrary, for parts that are more intricate, models using a retarding rate yield better results (RSC model and shear pARD-RSC model). In the sensor cover, it can be observed that different models yield the best results at different positions of the part. Close to the sprue, elongational flows are dominant, and the MRD and elongation models are the best. At Positions 1 and 2, the models using a retarding rate are superior. The only model that has a high prediction accuracy with an RMSD smaller than 0.1 in all points is the novel flow-dependent model since it is able to scale between flow regimes, as it applies a retarding rate in shear flows and neglects it in elongational flows.


**Table 8.** RMSD values between the experimental and simulated fiber orientation tensor in different geometries.
