*4.1. Simulation*

To set up a simulation that matches the SPR experiment, a number of fibers matching the experimental volume fraction was placed inside a shear cell. To guarantee a constant fiber volume fraction, Lees–Edwards periodic boundaries [42] were assigned to the cell faces perpendicular to the shearing direction (Figure 5). Mechanical properties were assigned to the fibers while rheological properties were assigned to the matrix. A simple shear flow field was applied to the matrix phase which translates into the hydrodynamic force term in the force balance calculation explained in Section 2.

**Figure 5.** Shear cell with periodic boundary conditions.

#### *4.2. Simulation Setup*

In order to obtain reliable rheological information from discontinuous fiber composites, accurate and repeatable initial conditions are needed [27]. To create the initial cluster of fibers for the simulation, the microstructure of the experimental samples was carefully characterized. The experimental values of orientation (aij) and fiber density (vol %) were reproduced by discretizing each through-thickness profile into a set number of layers as shown in Figure 6. For each layer the average orientation tensor components were used to reconstruct the fiber orientation state by using a Fourier Series expansion. It is worth noticing that as the thickness of the individual layers is decreased, to obtain better resolution, the target off-plane orientation becomes harder to achieve.

**Figure 6.** Through-thickness discretization for (**a**) fiber concentration and (**b**) fiber orientation of a compression molded sample. Experimental (black) and discretized, computational (red) data is shown [43].

The average FLs (LN, LW) were assigned as a global length distribution to the complete cluster of fibers. This distribution was obtained by fitting a Weibull probability distribution function to the experimental length measurement. The result was a heterogeneous cluster of fibers closely resembling the experimental condition (Figure 7). Prior to imposing the shear flow, a short simulation was conducted with no velocity field and only interaction forces present. This was applied as a relaxation step, so fibers are not overlapping or in extreme proximity at the start of the shearing simulation. This step does not modify the initial orientation state in a significant way.

**Figure 7.** Computational compression molded cluster, with x as the shearing direction (a11). Fibers are colored as a function of their XY (=12) planar orientation [43].

A simple shear flow with a rate of deformation of 1 s<sup>−</sup><sup>1</sup> was applied to the shear cell. The viscosity was set to 110 Pas using the properties of the neat PP. The dimension of the cell is dictated by the SPR gap in Z direction of 2 mm. To allow free rotation, the dimension in X and Y were set to 1.1 times the length of the longest fiber. The cluster properties corresponding to both sample preparation methods are listed in Table 2.


**Table 2.** Cluster properties for injection- and compression molded samples.

The simulation time was set to 60 s to obtain a total strain of 60. Cell walls parallel to the YZ plane have a periodic boundary condition. A tight array of fixed fibers is placed on the upper and lower boundaries to represent the SPR walls. The additional input parameters for the simulation are listed in Table 3. The simulation results were the coordinates of each fiber at every time step. With this information orientation tensors were calculated as a function of time and strain.


**Table 3.** Shear cell properties for injection- and compression molded samples.

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

#### *5.1. Injection Molded Samples*

IM provides a repeatable flow history that will yield to FO that is consistent among multiple samples [4,27]. The global diagonal components of the orientation tensor at zero shear strain are shown in Figure 8. The common core–shell structure associated with molded composites can be seen. The core layer consists of fibers predominantly aligned in cross-flow direction (a22), while fibers in the shell are oriented along the flow direction (a11) due to the fountain flow e ffect. The orientation in thickness direction a33 is uniform and generally low with average values less than 0.06 [4].

**Figure 8.** Initial fiber orientation for injection molded samples for simple shear flow tests.

The average orientation tensor as a function of shear strain is shown in Figure 9. For both simulation and experiment fibers gradually align in the flow direction over 60 shear strain units. The general trend of a11 for the experimental data begins at 0.53 and evolves to a value of 0.63. In contrast, the simulation starts at a value for a11 at 0.60 and evolves to a value of 0.74. The computational initial condition di ffers from the experimental one due to two main reasons: First, the complex core–shell structure present in the IM sample is di fficult to discretize (Figure 8). There will be a loss of accuracy by averaging the orientation data into seven layers. Second, the FL of the experimental sample is very high which increases the shear cell size and therefore computational time. Therefore, the FL was truncated in the simulation. As the initial conditions were not accurately reproduced, a di fferent rate of orientation evolution was obtained in the PLS. However, experiment and simulation still show similarities in their behavior. Both show a slight drop in a11 before fibers start aligning in shearing direction. This drop is less pronounced in the PLS, as fibers start o ff with a higher a11 alignment. An initial decrease in a11 was also found in [27]. In their work, this phenomenon was explained as follows. Fibers oriented slightly out of the 1–2 plane need to flip first to become oriented in flow direction, leading to an increase in a13, thus a decrease in a11. This increase in a13 could be seen in the experimental data. Depending on the initial value of a13, the alignment in flow direction is different. If a13 is positive, the fiber will continue to align in flow direction. However, an initial negative a13 orientation would require the fiber to flip before becoming aligned in flow direction. It is expected that a fiber with a negative a13 would take longer to align in flow direction [27]. A slight initial decrease in a11 was also reported by [44]. The used SPR was limited to 60 strain units and a steady state response associated with fiber suspensions could not be obtained.

**Figure 9.** Experimental (black) and predicted (red) fiber orientation evolution of injection molded samples.

The simulation results for a11 through thickness were smoothed and plotted for different strains in Figure 10. Under a constant shear, the complex heterogeneous orientation profile transitions into a more homogeneous profile.

**Figure 10.** Smoothed computational a11 evolution of injection molded samples at varying shear strains.

#### *5.2. Compression Molded Samples*

It has been shown that IM samples provide a repeatable initial FO; however, IM specimens show a complex core–shell structure, which is difficult to reproduce computationally (Figure 8). In addition, samples showed a high alignment in shearing direction at zero strain; therefore, only a low further alignment of a11 could be observed (Figure 9). CM samples with a controlled FL, FC, homogenous structure, and with a low a11 alignment were consequently manufactured to accurately reproduce the initial conditions (Figure 6).

The average orientation through sample thickness as a function of shear strain is shown in Figure 11. Initial conditions could be accurately reproduced computationally and fit the experimental data. The standard deviation for experimental values at zero strain is low, which demonstrates a repeatable initial FO, validating the CM sample preparation method.

**Figure 11.** Experimental (black) and predicted (red) fiber orientation evolution of compression molded samples [43].

The a11 component starts at 0.36 and transitions to a steady-state value of approximately 0.7 for both the experiment and the PLS. A steady state is reached at a total strain of 50. A slight initial drop in a11 could be seen in the experimental data, but was not present for the PLS. The simulation shows faster orientation evolution than the experiment. This phenomenon has also been reported in literature [45,46] for other diffusion models. It has been established that Jeffrey's hydrodynamic model and Folgar–Tucker's model predict a faster transient orientation rate when compared to related experiments [46]. Therefore, recent models have been developed with slow fiber orientation kinetics. Wang et al. [12] introduced the reduced strain closure model. This model employs a scalar factor κ to slow the orientation dynamics. Results still showed a quicker initial rise of flow-direction orientation, but achieved nearly the same steady state orientation as obtained with experiments [45]. In an attempt to slow the orientation kinetics predicted, Phelps et al. [13] derived the anisotropic rotary diffusion model.

The distortion of the flow field due to the presence of particles gives rise to the particular rheology of suspensions. As the aspect ratio (*rp*) of the particles increases, this effect becomes dependent on the orientation state, as described by Dinh et al. [47]. Lindstrom et al. [34] conducted a PLS with individual fibers in simple shear flow. They observed orbit periods were overestimated when not including two-way coupling. They concluded that particle–fluid interaction is essential to the fiber's dynamic behavior in shear flow. In that regard, neglecting the two-way coupling between the fiber orientation state and the underlying flow field in a PLS, can be a considerable source of error. However, various authors have shown that including two-way coupling has marginal effect on the orientation state for flow in common geometries. Mezi et al. [48] pointed out that accounting for the coupling effect had little impact on the fiber orientation distribution for the fluid flow in a planar channel, while having a considerable effect in the pressure drop. They also observed coupling increased the size of corner vortex in a contraction geometry. This last effect was also noted by VerWeyst et al. [49], who observed

moderate impact of coupling on the fiber orientation in a center-gated disk. They showed effects of coupling decayed rapidly with increasing distance from the center of the disk. The difficulty with these approaches lies in the large computational domain required to resolve both the spatial and orientation domains. With PLS the computational cost increases rapidly as the two-way coupling must be calculated for each fiber.

In a PLS study for SFTs under simple shear, Kugler et al. [50] also observed a faster alignment in the simulation with respect to the experimental data. They suggested a reduced shear strain to correct for the lack of two-way coupling. Similarly, in an experimental study of fiber attrition employing a Couette rheometer, Moritzer et al. [51] use the scale factor (1 − φ)) to account for the reduced rate of deformations experienced by fibers in the suspension.

Published work indicates that the coupling effect might have a greater impact on the fiber orientation in the dynamic portion of the orientation evolution, and that aspect ratio and volume fraction are scaling factors of such effect. A PLS study to determine quantitatively the impact of coupling under simple shear flow remains to be done; one that compares the predicted orientation to the complete microstructural data determined experimentally.

The a13 and a12 components of orientation are shown in Figure 12. In both prediction and experiment a13 and a12 show close agreemen<sup>t</sup> with the final steady state value. The measured a13 component shows small oscillations. This behavior indicates some fibers orbit in the plane of shear. The predicted a13 component also shows these orbits, however, with a shorter period. The initial measured a13 component has a negative value of −6 × <sup>10</sup>−4, while the computational cluster started with a very small, ye<sup>t</sup> positive, value of 6 × 10−5. Literature suggests that the initial value of a13 impacts whether components of orientation evolve monotonically to a steady state or present an overshoot and/or undershoot prior to reaching a steady state [52]. The fact that the initial computational cluster differed from the experimental one in this aspect, is a factor in why the a11 component oriented directly to steady state in the simulation, while the experimental a11 did show a slight undershoot (Figure 11). It can be seen that a12, for the simulation, increases continuously with shear strain. As fibers oriented in the 1–2 plane orient themselves in 1, the value of a12 should approach zero.

**Figure 12.** Experimental (black) and predicted (red) off-diagonal orientation components (**a**) a12 and (**b**) a13 as a function of shear strain.

The mentioned faster orientation kinetics in the PLS are also shown in Figure 13. Compared to the initial orientation, the experimental core–shell structure is largely unchanged at 12.5 strain units. At the same applied strain, the simulation already shows a large transition to a rather homogenous structure. At 60 strain units, the core–shell structure disappears, a homogeneous profile is reached and a good match between experiment and the PLS was observed. It is assumed that shearing beyond 60 strain units would lead to a steady state FO through thickness [27]. The smoothed computational

a11 evolution as a function of shear strain is shown Figure 14. With an increase in shear strain, the heterogeneous orientation profile transitions into a homogeneous steady state. The most significant change in orientation occurs in the core of the sample. The rate of change slows down as the orientation approaches steady state.

**Figure 13.** Experimental (black) and predicted (red) a11 values averaged through sample thickness at varying shear strains for compression molded samples. (**a**) Initial fiber orientation distribution at 0 shear strain, (**b**) fiber orientation distribution after 5 strain units, (**c**) 12.5, (**d**) 20, (**e**) 40, and (**f**) 60 shear strain units [43].

**Figure 14.** Smoothed computational a11 evolution of compression molded samples [43].

With the proposed CM technique, it was possible to manufacture samples with low alignment in shearing direction, which allowed a larger orientation transition. Higher alignment and faster orientation evolution was observed for the CM samples when compared to the IM ones. The faster alignment can be attributed to the lower FL of the CM samples. The long fibers present in the IM samples require more energy to rotate and also hinder the motion of adjacent fibers.

#### **6. Conclusions and Outlook**

A PLS was used to simulate the FO evolution of glass fiber-reinforced PP plates sheared in a SPR. The PLS results showed a faster orientation evolution at the beginning of the shearing process compared to the experimental data. However, an agreemen<sup>t</sup> with the final orientation state was observed for the CM plates after fine-tuning the initial condition discretization method. In the experiment, an early decrease in a11 was observed for IM and CM samples. The same behavior has been reported in literature by [27]. However, the author used shorter fibers compared to the FL in this work. It appears that this phenomenon is more easily observed with longer fibers due to a slower rate of FO.

A reliable sample preparation method for suspension rheology was developed. Repeatable and controlled initial orientation can be achieved through the presented CM technique. Even though it has been shown in [27] that CM samples do not exhibit repeatable orientation evolution data during the startup of simple shear flow in a SPR, the CM technique employed in this work di ffers significantly by manually depositing individual and aligned strands.

As the dynamics of fiber suspensions are not ye<sup>t</sup> well captured by the PLS, future work will focus on implementing the e ffect of fiber fluid coupling and studying its impact under di fferent combinations of the dimensionless number φ*rp*.

**Author Contributions:** S.A.S. and A.B.S. conceived and designed the experiments. S.A.S. was responsible for performing the experimental studies and the microstructure analysis. A.B.S. performed the simulations. S.A.S. and A.B.S. analyzed the data and wrote the paper. T.O. supervised the project and was involved in all stages of the research. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was supported by the National Science Foundation under Grant No. 1633967.

**Acknowledgments:** The authors would like to thank SABIC ® for their ongoing support and collaboration with our research group and their material supply. Additionally, we would also like to thank our colleagues from the Polymer Engineering Center who provided insight and expertise that greatly assisted the research.

**Conflicts of Interest:** The authors declare no conflicts of interest.
