Experimental Validation of a Direct Fiber Model for Orientation Prediction

: Predicting the ﬁber orientation of reinforced molded components is required to improve their performance and safety. Continuum-based models for ﬁber orientation are computationally very e ﬃ cient; however, they lack in a linked theory between ﬁber attrition, ﬁber–matrix separation and ﬁber alignment. This work, therefore, employs a particle level simulation which was used to simulate the ﬁber orientation evolution within a sliding plate rheometer. In the model, each ﬁber is accounted for and represented as a chain of linked rigid segments. Fibers experience hydrodynamic forces, elastic forces, and interaction forces. To validate this fundamental modeling approach, injection and compression molded reinforced polypropylene samples were subjected to a simple shear ﬂow using a sliding plate rheometer. Microcomputed tomography was used to measure the orientation tensor up to 60 shear strain units. The fully characterized microstructure at zero shear strain was used to reproduce the initial conditions in the particle level simulation. Fibers were placed in a periodic boundary cell, and an idealized simple shear ﬂow ﬁeld was applied. The model showed a faster orientation evolution at the start of the shearing process. However, agreement with the steady-state aligned orientation for compression molded samples was found.


Introduction
Over the past decade, the use of long fiber-reinforced thermoplastics (LFTs) has gained wide acceptance in the automotive industry in order to meet the increasingly tightened corporate average fuel efficiency standards [1,2]. As more efficient engines and electric powertrains will not carry the whole load, the automotive industry is pressured into finding alternative solutions. One of those solutions has been light weighting car components by retrofitting existing structures made out of steel [3][4][5][6][7]. LFTs are being considered as a substitute because they show high performance in terms of mechanical properties, are lightweight, noncorrosive, and can be tailored to satisfy performance requirements [3,5,6,8,9].
Although the mechanical advantages of LFTs are widely reported in literature, the uncertainty remains being able to accurately model this heterogenous material class. It has been shown by numerous researchers that LFTs' mechanical performance and dimensional stability are a direct function of their microstructure, and therefore depend on fiber orientation (FO), fiber length (FL), and fiber concentration (FC) [7]. Parts are stronger in the direction of the fiber alignment and if fiber length and volume fraction (ϕ) are increased [10]. Being able to accurately predict the microstructure of molded components is a key factor in the automotive industry, not only for design calculations and a successful process, but also for an optimized part's performance and for guaranteeing parts are safely introduced into vehicles. In the model, inertial effects are neglected due to the low Reynolds numbers (Re) that result from the high viscosity of the polymer matrix. It has been shown experimentally by Hoffman [28] and Barnes [29] that inertial effects become important at Re ≥ 10 −3 . In our work, the Re was calculated to be in a range of 10 −9 -10 −4 . Long-range hydrodynamic interactions are neglected as well due to the fluid's high viscosity [30]. Additionally, extensional and torsional deformations are neglected as fluid forces are not sufficiently high to cause fiber stretching or torsional deformation. Only bending deformation is included in the model. Brownian motion can be neglected as the Peclet number (Pe) is well above 10 3 and at this point Brownian interactions are disrupted by hydrodynamic forces [31]. Buoyant effects are neglected in the model as well [7,32]. The presence of fillers in any suspension modifies the flow field, especially in the semi-diluted and more concentrated regimes. However, the coupling between particle and fluid is not considered in the present formulation due to its high computational cost [32].
The force calculation for a segment includes the drag force , the interaction force with an adjacent segment , and intra-fiber forces exerted by internal connection loads at the nodes. The translation equation of motion is The rotational equation of motion is equally derived and includes elastic recovery terms and a hydrodynamic torque where describes the shortest distance vector between two segments. Some of these quantities are presented in Figure 2. A fiber divided into more than one element requires an extra constraint that enforces connectivity between the different elements 0 = + × ( +1 − ) − +1 (3) Figure 1. Particle-level simulation: modeling single fiber and macroscale interactions [7].
In the model, inertial effects are neglected due to the low Reynolds numbers (Re) that result from the high viscosity of the polymer matrix. It has been shown experimentally by Hoffman [28] and Barnes [29] that inertial effects become important at Re ≥ 10 −3 . In our work, the Re was calculated to be in a range of 10 −9 -10 −4 . Long-range hydrodynamic interactions are neglected as well due to the fluid's high viscosity [30]. Additionally, extensional and torsional deformations are neglected as fluid forces are not sufficiently high to cause fiber stretching or torsional deformation. Only bending deformation is included in the model. Brownian motion can be neglected as the Peclet number (Pe) is well above 10 3 and at this point Brownian interactions are disrupted by hydrodynamic forces [31]. Buoyant effects are neglected in the model as well [7,32]. The presence of fillers in any suspension modifies the flow field, especially in the semi-diluted and more concentrated regimes. However, the coupling between particle and fluid is not considered in the present formulation due to its high computational cost [32].
The force calculation for a segment i includes the drag force F H i , the interaction force with an adjacent segment j F C ij , and intra-fiber forces X i exerted by internal connection loads at the nodes. The translation equation of motion is The rotational equation of motion is equally derived and includes elastic recovery terms M b and a hydrodynamic torque where r ij describes the shortest distance vector between two segments. Some of these quantities are presented in Figure 2. A fiber divided into more than one element requires an extra constraint that enforces connectivity between the different elements Hydrodynamic effects are approximated by modeling each fiber segment as a chain of beads. The hydrodynamic force F H i is the sum of forces experienced by the beads F H k and is described as F H k is given by Stokes law, where k describes the number of beads, a the bead radius, µ the polymer matrix viscosity, U ∞ k the surrounding fluid velocity, and u k the velocity of the bead k (u k = u i + ω i × r k ).
F H i is then be written as The rotational equation of motion is equally derived and includes elastic recovery terms and a hydrodynamic torque where describes the shortest distance vector between two segments. Some of these quantities are presented in Figure 2. A fiber divided into more than one element requires an extra constraint that enforces connectivity between the different elements Due to the fluid's vorticity each bead is also subjected to a hydrodynamic torque T H where the hydrodynamic contribution T H k of bead k is Ω ∞ k is the vorticity of the surrounding fluid and ω k is the angular velocity of the bead k. Substituting the expression of T H k and the expression of u k into the fluid's vorticity we can write The fiber-fiber interaction force F C ij is decomposed in a normal force F N ij and a tangential force F C ij describes the friction between the segments and starts acting when the distance between adjacent fibers, d ij , is below a defined threshold d threshold . F N ij is an excluded volume force which is implemented as a discrete penalty method [7,32,33]. F N ij usually has a function of the following type [26,34,35], where d ij is the shortest distance between rod i and rod j, D the fiber diameter, and n ij is the vector along the closest distance between the rods. A and B are parameters for which different values have been presented in the literature [26,34]. In this work, A has been chosen empirically to avoid the overlapping of fibers. A value of 1000 N was the minimum constant value at which penetrations cease to occur in the simulation [32]. B was set to 2. The friction force F T ij between segment i and j is computed as µ f is the Coulomb coefficient between fibers, ∆u ij is the vector of the relative velocity between elements.
Bending of a fiber was approximated by using elastic beam theory with M b i as the bending moment, α as the angle between segments, and l as the segment length. A linear system of equations was assembled and solved to find the velocities and connective forces in each fiber. These velocities were integrated over time using an explicit Euler scheme to determine the fiber trajectory during the simulation [7].

Sample Preparation
The fiber suspension used in this work was a 10 and 20 wt % glass fiber in a polypropylene (PP) matrix. The 20 wt % (STAMAX PPGF20) was commercially available and provided by SABIC TM . The material was supplied in the form of coated pellets with a nominal length of 15 mm, which also represents the initial and uniform length of the glass fibers. The used E-glass fibers ( = 2.55 g/cm 3 ) are chemically coupled to the PP matrix ( = 0.91 g/cm 3 ). The fiber diameter was measured to be 19 ± 1 µm using an optical microscope. The 10 wt % was achieved by mixing higher FCs with neat PP (SABIC TM PP579S) in a cement mixer. The neat PP is identic to the matrix material for the commercial available STAMAX pellets. The matrix was considered as a generalized Newtonian fluid under the test conditions presented. The authors of [27,36] showed that the presence of polymer matrix eliminates any concerns regarding fiber sedimentation due to gravity.
Two separate processes were used to prepare samples that were tested in the SPR. The first method combined extrusion to disperse fibers and compression molding to form proper sample dimensions. The second approach used IM to produce plates.
The first sample preparation method used an Extrudex Kunstoffmaschinen single-screw extruder to process the pultruded material. The 45 mm 30 L/D extruder was equipped with a gradually tapering screw and fitted with a 3 mm die. The temperature zones of the extruder were set to 210, 210, 220, 220, 230, 230, and 230 • C and the die temperature was set to 230 • C. The composite was extruded at 5 rpm. Due to the low processing speed, most of the initial FL was maintained in the extruded strand (L N = 6.7 mm, L W = 11.4 mm). Computational time for the mechanistic model simulation is geometrically proportional to the maximum detected FL. To reduce computation, the strands were pelletized to 3.2 mm. Initial compression molding trials showed that manual alignment of pellets in the mold did not yield a repeatable initial FO, because pellets could move easily and rotate. Therefore, pellets were re-extruded and strands were cut to fit mold dimensions. The mold geometry was a rectangular prism (14 mm × 14 mm × 2.1 mm). Analysis of the re-extruded strands showed a L N of 0.83 mm and L W of 1.53 mm. The mold was coated and placed on a heating plate to cause partial melting of the aligned strands. This facilitated space reduction between irregular shaped strands. The plate was extracted, flipped 180 • , and re-molten to remove trapped air bubbles between strands. Failure to remove caught air bubbles would alter matrix rheology. Plates were compression molded at 210 • C and a load of 1000 lbs. Final fiber properties are summarized in Table 1.  [-] 0.86 0.60 a 22 [-] 0.11 0.37 a 33 [-] 0.03 0.03 The second sample preparation procedure involved IM a simple plate geometry (102 × 305 × 2.85 mm 3 ). The cavity was filled through a 20 mm edge-gate with the same thickness as the plate and was fed through a 17 mm full-round runner. Parts were molded on a 130 ton Supermac Machinery SM-130 IM machine. The melt temperature was set to 250 • C, and back and holding pressure to 5 and 300 bar, respectively. Injection and holding time were set to 2 and 22 s, respectively. A full microstructure analysis was conducted. Fiber properties are summarized in Table 1. Samples for the SPR were only extracted from regions which showed a clear developed shell-core structure and identic FCs [4].
For both sample preparation methods, the SPR sample sizes were calculated by adding 20 mm to the stroke applied to be able to extract pure shear samples for FO analysis.

Sliding Plate Rheometer
FO evolution data was obtained in a SPR by shearing samples in a controlled simple shear flow ( Figure 3). The rheometer was based on the design of [37]. The SPR was contained in a convection oven and the rectilinear plate displacement was generated by an Interlaken 3300 universal testing instrument. The rheometer had an effective surface of 100 × 300 mm 2 and a maximum stroke of 120 mm. The SPR gap thickness could be altered. However, it was chosen to be 2 mm, as this work yields to aid in the understanding of how fibers flow and orient themselves during the IM process. Typically, IM parts show a thickness of 2 mm. The maximum possible displacement and the chosen gap size, limited the deformation that could be imposed on the sample to a shear strain of 60. Both, shear rate and shear strain, were programmable through the Wintest ® Software (Bose Corporation, Eden Prairie, MN, USA).  Table 1. Samples for the SPR were only extracted from regions which showed a clear developed shell-core structure and identic FCs [4]. For both sample preparation methods, the SPR sample sizes were calculated by adding 20 mm to the stroke applied to be able to extract pure shear samples for FO analysis.

Sliding Plate Rheometer
FO evolution data was obtained in a SPR by shearing samples in a controlled simple shear flow ( Figure 3). The rheometer was based on the design of [37]. The SPR was contained in a convection oven and the rectilinear plate displacement was generated by an Interlaken 3300 universal testing instrument. The rheometer had an effective surface of 100 × 300 mm 2 and a maximum stroke of 120 mm. The SPR gap thickness could be altered. However, it was chosen to be 2 mm, as this work yields to aid in the understanding of how fibers flow and orient themselves during the IM process. Typically, IM parts show a thickness of 2 mm. The maximum possible displacement and the chosen gap size, limited the deformation that could be imposed on the sample to a shear strain of 60. Both, shear rate and shear strain, were programmable through the Wintest ® Software (Bose Corporation, Eden Prairie, MN, USA). The experimental procedure was based on the SPR experiment conducted by [27]. The rheometer was heated to 260 °C for 2 h prior to sample loading. Upon loading, the test specimen was rotated 90° with respect to the extrusion/filling direction, effectively swapping the a11 and the a22 components of the orientation tensor. This allowed for a low alignment in the shearing direction, so a larger change in orientation could be observed. The sample was secured between the rheometer plates and allowed to melt evenly before the plates were tightened to a final gap of 2 mm (=initial condition). As the initial thickness of the test specimen was larger than the final gap thickness, samples were slightly compressed when tightening the plates to guarantee full contact. After an additional 10 min of heating, the sample was sheared at a rate of 1 s −1 . Forced convection was used to cool the sample to preserve its shape and FO for further analysis. Five repetitions per testing condition were used to ensure accuracy and repeatability of results.

Measurement of Fiber Microstructure
A fully characterized microstructure was required to accurately reproduce the initial conditions 2 1 3 Figure 3. Sliding plate rheometer and defined coordinate system.
The experimental procedure was based on the SPR experiment conducted by [27]. The rheometer was heated to 260 • C for 2 h prior to sample loading. Upon loading, the test specimen was rotated 90 • with respect to the extrusion/filling direction, effectively swapping the a 11 and the a 22 components of the orientation tensor. This allowed for a low alignment in the shearing direction, so a larger change in orientation could be observed. The sample was secured between the rheometer plates and allowed to melt evenly before the plates were tightened to a final gap of 2 mm (=initial condition). As the initial thickness of the test specimen was larger than the final gap thickness, samples were slightly compressed when tightening the plates to guarantee full contact. After an additional 10 min of heating, the sample was sheared at a rate of 1 s −1 . Forced convection was used to cool the sample to preserve its shape and FO for further analysis. Five repetitions per testing condition were used to ensure accuracy and repeatability of results.

Measurement of Fiber Microstructure
A fully characterized microstructure was required to accurately reproduce the initial conditions in the direct fiber simulation ( Figure 4). FO was determined by using the X-ray microcomputed tomography (µCT) approach. A sample of dimensions 10 × 10 × 2 mm 3 was loaded into a sample holder and then placed on a rotating platform. The X-rays penetrated the sample and were absorbed differently depending on the configuration of the sample's constituents. The detector recorded the attenuated X-rays as radiographs at incremental angles during the rotation of the sample to achieve a full scan of the sample. After completion of the 3D reconstruction, the µCT data set was processed with the VG StudioMAX (Volume Graphics) software to quantify the FO distribution using the structure tensor approach [4,38,39]. Samples were scanned with an industrial µCT system (Metrotom 800, Carl Zeiss AG, Oberkochen, Germany). Throughout the reported studies, the voltage was set to 80 V, the current to 100 A, the integration time to 1000 ms, the gain to 8, the number of projections to 2200, and the voxel size was set to 5 µm.
J. Compos. Sci. 2020, 4, x 7 of 17 800, Carl Zeiss AG, Oberkochen, Germany). Throughout the reported studies, the voltage was set to 80 V, the current to 100 A, the integration time to 1000 ms, the gain to 8, the number of projections to 2200, and the voxel size was set to 5 μm. FC was quantified with VG StudioMAX as well. The μCT data set was converted into a stack of 2D cross-sectional images aligned normal to the thickness direction. The grayscale images were transformed into binary images by thresholding, which separated the image into black (matrix) and white (fibers) pixels. Subsequently, the fiber volume fraction through the thickness was calculated [4].
The FL measurement technique presented in [40] was employed in this work. This technique consists of fiber dispersion and a fully automated image processing algorithm to quantify the fiber length distribution (FLD). It has been shown in [41] that downsampling methods yield to a preferentially capture of long fibers, and thus skew the real FLD. In this work, a correction was therefore applied to all results as published by [41].

Simulation
To set up a simulation that matches the SPR experiment, a number of fibers matching the FC was quantified with VG StudioMAX as well. The µCT data set was converted into a stack of 2D cross-sectional images aligned normal to the thickness direction. The grayscale images were transformed into binary images by thresholding, which separated the image into black (matrix) and white (fibers) pixels. Subsequently, the fiber volume fraction through the thickness was calculated [4].
The FL measurement technique presented in [40] was employed in this work. This technique consists of fiber dispersion and a fully automated image processing algorithm to quantify the fiber length distribution (FLD). It has been shown in [41] that downsampling methods yield to a preferentially capture of long fibers, and thus skew the real FLD. In this work, a correction was therefore applied to all results as published by [41].

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.

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.
(a) (b) 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

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 (a ij ) 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.

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.
(a) (b) 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 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 (L N , L W ) 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. A simple shear flow with a rate of deformation of 1 s −1 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. 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.  A simple shear flow with a rate of deformation of 1 s −1 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. 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.

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 (a 22 ), while fibers in the shell are oriented along the flow direction (a 11 ) due to the fountain flow effect. The orientation in thickness direction a 33 is uniform and generally low with average values less than 0.06 [4].

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 effect. The orientation in thickness direction a33 is uniform and generally low with average values less than 0.06 [4]. 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 differs from the experimental one due to two main reasons: First, the complex core-shell structure present in the IM sample is difficult 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 different 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 off 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 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 a 11 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 a 11 at 0.60 and evolves to a value of 0.74. The computational initial condition differs from the experimental one due to two main reasons: First, the complex core-shell structure present in the IM sample is difficult 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 different 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 a 11 before fibers start aligning in shearing direction. This drop is less pronounced in the PLS, as fibers start off with a higher a 11 alignment. An initial decrease in a 11 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 a 13 , thus a decrease in a 11 . This increase in a 13 could be seen in the experimental data. Depending on the initial value of a 13 , the alignment in flow direction is different. If a 13 is positive, the fiber will continue to align in flow direction. However, an initial negative a 13 orientation would require the fiber to flip before becoming aligned in flow direction. It is expected that a fiber with a negative a 13 would take longer to align in flow direction [27]. A slight initial decrease in a 11 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. 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.

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

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,

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 a 11 could be observed (Figure 9). CM samples with a controlled FL, FC, homogenous structure, and with a low a 11 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. 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 ( ) 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 Figure 11. Experimental (black) and predicted (red) fiber orientation evolution of compression molded samples [43].
The a 11 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 a 11 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 (r p ) 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 a 13 and a 12 components of orientation are shown in Figure 12. In both prediction and experiment a 13 and a 12 show close agreement with the final steady state value. The measured a 13 component shows small oscillations. This behavior indicates some fibers orbit in the plane of shear. The predicted a 13 component also shows these orbits, however, with a shorter period. The initial measured a 13 component has a negative value of −6 × 10 −4 , while the computational cluster started with a very small, yet positive, value of 6 × 10 −5 . Literature suggests that the initial value of a 13 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 a 11 component oriented directly to steady state in the simulation, while the experimental a 11 did show a slight undershoot ( Figure 11). It can be seen that a 12 , for the simulation, increases continuously with shear strain. As fibers oriented in the 1-2 plane orient themselves in 1, the value of a 12 should approach zero. J. Compos. Sci. 2020, 4, x 13 of 17 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 agreement 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 × 10 −4 , while the computational cluster started with a very small, yet 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. 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 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 a 11 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.  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  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 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.

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 agreement 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 a 11 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 differs significantly by manually depositing individual and aligned strands.
As the dynamics of fiber suspensions are not yet well captured by the PLS, future work will focus on implementing the effect of fiber fluid coupling and studying its impact under different combinations of the dimensionless number φr p .