A Flow-Dependent Fiber Orientation Model

: The mechanical performance of ﬁber reinforced polymers is dependent on the process-induced ﬁber orientation. In this work, we focus on the prediction of the ﬁber orientation in an injection-molded short ﬁber reinforced thermoplastic part using an original multi-scale modeling approach. A particle-based model developed for shear ﬂows is extended to elongational ﬂows. This mechanistic model for elongational ﬂows is validated using an experiment, which was conducted for a long ﬁber reinforced polymer. The inﬂuence of several ﬁber descriptors and ﬂuid viscosity on ﬁber orientation under elongational ﬂow is studied at the micro-scale. Based on this sensitivity analysis, a common parameter set for a continuum-based ﬁber orientation macroscopic model is deﬁned under elongational ﬂow. We then develop a novel ﬂow-dependent macroscopic ﬁber orientation, which takes into consideration the effect of both elongational and shear ﬂow on the ﬁber orientation evolution during the ﬁlling of a mold cavity. The model is objective and shows better performance in comparison to state-of-the-art ﬁber orientation models when compared to µ CT-based ﬁber orientation measurements for several industrial parts. The model is implemented using the simulation software Autodesk Moldﬂow Insight Scandium R (cid:13) 2019.


Introduction
The current needs for weight reduction and eco-friendly products (e.g., natural fiber composites) have fueled the use of polymers in various industries. They are of particular importance in the automotive and aerospace industries due to their low weight and durability. The polymers in most cases are further reinforced with fibers to enhance their mechanical properties. Those are strongly influenced by fiber volume fraction, fiber length distribution, and fiber orientation.
Fiber reinforced polymers show typically anisotropic material properties. The local orientation of fibers depends on the manufacturing process and the complexity of the part being produced. The prediction of the failure and lifetime of a part using simulative methods instead of experimental methods can significantly help to reduce the costs and speed up the development of new parts. The final fiber orientation induced by the manufacturing process is an important parameter for the structural simulation of the part. The current state-of-the-art fiber orientation models based on Jeffrey's equation [2] have been developed and implemented in different commercial software. These seek to emulate several effects like strain reduction [3][4][5] and anisotropic rotary diffusion [6][7][8] on fiber orientation. All these models require phenomenological parameters derived from experimental or microsimulation results to obtain a suitable prediction accuracy.
The experimental determination of the phenomenological parameter to describe the steady state fiber orientation by diffusion was for example conducted by Folgar and Tucker or Bay and Tucker [9,10]. In order to determine phenomenological parameters, which describe the fiber orientation evolution by the retarding rate, it is necessary to measure the transient fiber orientation development. Consequently, the methods are more demanding. Stover et al. [11] measured transient fiber orientations in semi-dilute solutions using a Couette device, a tracer fiber, and video cameras. For highly concentrated solutions, in situ methods with optical methods are not applicable yet. To overcome this shortcoming, experiments in homogeneous flows with a repeatable starting condition have been developed. Eberle et al. [12] measured fiber orientation evolution for a 30% wt. short fiber reinforced polybutylene terephthalate (PBT) in a cone and plate rheometer with specially shaped donut-like samples. Ortman et al. [13] used a sliding plate rheometer, where it is possible to control the initial conditions [14], to measure fiber orientation evolution. Kugler et al. [15] used the same setup to determine phenomenological parameters for a short fiber reinforced PBT-GF30. Recently, Perumal et al. [16] measured transient fiber orientation evolution for a 30% wt. glass fiber reinforced Nylon-6 in a parallel plate rheometer. All stated approaches determine fiber orientation parameters in shear flow. Lambert et al. [17,18] determined fiber orientation parameters for the first time in elongational flow.
Additionally to experimental determination, microscopic fiber simulation can be used to evaluate macroscopic fiber orientation parameters. Modeling approaches on the microscopic scale have the aim to approximate the physical behavior of the composite more accurately. Many authors determined the parameters based on microscopic simulation, for example Mezher et al. [19] and Perez [1]. For a detailed review, refer for example to [20].
In this work, we use a multi-scale simulation chain to enhance the fiber orientation prediction for complex industrial parts. The simulation workflow consists of three steps. Firstly, a virtual flow test on the particle scale (micro) is conducted. Secondly, macroscopic fiber orientation model parameters are determined based on the resulting fiber orientation evolution from the virtual flow test. Thirdly, the final fiber orientation in a part is obtained using continuum-based macro models with the material-dependent optimal parameters. The multi-scale process simulation is sketched in Figure 1. Current fiber orientation models, on both scales, are mostly focused on the effects of shear flow. However, experiments in elongation flow [17,18] showed that fiber orientation develops differently under elongational flows. On the other hand, during injection molding, complex flow fields combining shear and elongational flow usually take place. Considering only shear-fitted fiber orientation models may lead to discrepancies in the final fiber orientation predictions in a part. Recently, Chen et al. [21,22] introduced a flow-dependent strain reduction factor for the different state-of-the-art macroscopic fiber orientation models. This last model approach is able to differentiate the orientation evolution speed between shear, elongational, and rotational flow.
In this work, we focus primarily on the influence of elongational flow on the fiber orientation phenomenon using a simulation at the particle level. Finally, a novel flow-dependent fiber orientation model, scaling between shear and elongational flows, is proposed and implemented in Autodesk Moldflow Insight Scandium R 2019. In comparison to the approach proposed by Chen et al. [21,22], we propose a more general model since it not only considers the change in strain reduction, but also in diffusion.

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.

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. The components of p can be represented in terms of θ and φ as: A probability distribution function (PDF) ψ(p, t) can be used to represent the orientation of a population of fibers in space and time [9]. The PDF ψ(p, t)dp gives the probability of a fiber being directed between p and p+dp at a time t. The properties of the PDF are: 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 = p i p j ψ(p)dp (8) a ijkl = p i p j p k p l ψ(p)dp (9) a i...n = p i . . . p n ψ(p)dp.
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 = a ij and A = a ijkl 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], The fourth order tensor is symmetric for any pair of indices: and it includes the entire information contained in the second order tensor [23]: 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].

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.
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+λ 2 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]: whereȦ h corresponds to the hydrodynamic effects andȦ d is the diffusion term, which takes care of the fiber-fiber interactions. The term C I is the interaction coefficient, andγ = √ 2D : 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][25][26], exact closures [27][28][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: 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: 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 C I 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.
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: 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: 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],

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 segment 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 ∞ . 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 segment i are: hydrodynamic forces F H i from the fluid, the interaction force with neighboring segment j F C ij , and intra-fiber force from the adjacent segment i X i . M b i and M b i+1 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).
where r i is the segment vector and d ij 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. 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)).
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]: where N is the total number of fibers in the system, t is time, and p n (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.

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.

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: 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 v upper w and v lower w towards the negative y-direction and the positive y-direction, respectively. They are given according to: 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 v w|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 −4 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. 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.

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].

Material Sample Thickness Temperature Strain Rate
Max Strain Total Time 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 2 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. 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 agreement 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.  [17] in points and simulation data in continuous lines.

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 mm 2 . The simulation parameters are listed in Table 4. The RVE was compressed in the y-direction at a constant strain rate during 40 s. 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. The good agreement 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.

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 −1 and 100.0 s −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. 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 −4 m, 8 × 10 −4 m, and 1.5 × 10 −3 m. It can be observed from Figure 8 that the fiber length had a negligible effect on the fiber orientation. 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 −4 m and a radius of 5 × 10 −6 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 2 and a height of 3 × 10 −3 m. The volume fraction of each RVE is listed in Table 5. An elongational rate of 0.05 s −1 was applied until reaching a strain value of two. The fiber orientation evolution over strain is depicted in Figure 9.  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 −1 , a volume fraction of 18%, a fiber radius of 5 × 10 −6 m, and a fiber length of 2.5 × 10 −4 m were used. The RVE height was 3 × 10 −3 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. 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.

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. 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 C I = 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.

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 R 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.

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 M z [35]. Instead of the Manas-Zloczower number, the objective measure β by Astarita [36] can also be used.
The Manas-Zloczower number is defined by: 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 M z is objective. The second invariant of the rate of strain tensor is defined by tr D 2 [37], so: The second invariant of the vorticity tensor is defined as tr W 2 = 1 2 Ω 2 , where Ω = ∇u is the vorticity vector [37]. Therefore, the magnitude of the vorticity tensor is defined by: The Manas-Zloczower number characterizes the flow type in the following way: 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 z . If at a given point during the flow, M z < 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:Ȧ In pure shear, this equalsȦ =Ȧ s , and in pure elongational flow,Ȧ =Ȧ e . The derivativesȦ s andȦ 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 C I has to be fitted.
This leads to the following definition of the flow-dependent fiber orientation model: This linearly interpolated fiber orientation model was implemented in Autodesk Moldflow Insight Scandium R 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.

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.

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. Table 7. Parameters of the flow-dependent pARD-RSC model for a PBT -GF30.

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: • pARD-RSC model with parameters from shear flow only • pARD-RSC model with parameters from elongational flow only • Flow-dependent pARD-RSC model.
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 R 2019. The RSC and MRD models are pre-defined in Autodesk Moldflow R . 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 M z after 0.41 s of filling. For information, the cavity was completely filled after 0.63 s. We could observe that M z exhibited higher values at the core of the sprue and plate. The value of M z 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. 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.

Conclusions
We propose a user-friendly multi-scale virtual workflow for estimating the fiber orientation of an injection-molded fiber reinforced thermoplastic part. To do this, we use a particle-based mechanistic model, which is able to evaluate the fiber orientation under shear flow and elongational flow. We validate the mechanistic model under elongational flow in front of experimental data using long fibers [17,18]. The estimations of the mechanistic model are in agreement with the experimental results. By exploiting the mechanistic model, we find that the fiber orientation evolution under elongational flow is independent of the fiber length, the fiber volume content, and the elongational rate for short fibers. There is however a slight sensitivity to variations of the matrix polymer viscosity. In future work, it would be interesting to study the influence of fiber flexibility on the fiber orientation phenomenon, in particular for long fibers.
At the macro-scale, we introduce a novel flow-dependent fiber orientation model that is able to adjust the fiber orientation evolution as a function of the local flow type by using an objective scaling parameter: the Manas-Zloczower number. This is reached by defining flow specific retarding and anisotropic diffusion parameters. We implement the 3D flow-dependent fiber orientation model in Moldflow R using a Solver API feature. The fiber orientation results from our flow-dependent fiber orientation model are compared to µCT scans in three different injection-molded parts. The novel model gives not only comparable, but in some cases, more accurate results than current existing models. Additionally, it is the only model, that provides good results in a plate and a complex part. This shows that the flow dependency provides a more general modeling approach in front of the existing macroscopic fiber orientation models. We believe that the global performance of the model can still increase by incorporating fiber orientation-viscosity coupling in the macroscopic simulation. In our flow-dependent model, the fiber orientation evolution speed scales linearly between shear and elongation flows. A future enhancement could be the use of a different scaling approach and the incorporation of rotational flows. This novel multi-scale simulation approach consisting of the mechanistic model and the flow-dependent macro-model is significantly faster when compared to workflows with experimental parameter identification.

Patents
This work is partially covered by Patent No. 102020204586.0.