Experimental Validation of a Direct Fiber Model for Orientation Prediction

Abstract

Predicting the fiber orientation of reinforced molded components is required to improve their performance and safety. Continuum-based models for fiber orientation are computationally very efficient; however, they lack in a linked theory between fiber attrition, fiber–matrix separation and fiber alignment. This work, therefore, employs a particle level simulation which was used to simulate the fiber orientation evolution within a sliding plate rheometer. In the model, each fiber 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 flow 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 flow field 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.

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

Continuum-based models for fiber orientation [11,12,13,14], length distribution [15,16], and fiber content prediction [17,18,19] have been developed in the last decades. These models employ a probabilistic approach to obtain the fiber configuration evolution during processing which is based on the input parameters, such as fiber characteristics and flow conditions. These models are computationally very efficient and have been implemented into commercial software. However, they can only describe the mechanism of fiber–matrix separation, fiber-attrition, and fiber alignment individually. They lack a linked theory between the three mechanisms and their interdependency [7]. Additionally, these models require experimentally determined fitting parameters. Conducting these experiments is costly, lengthy, and their provided information is limited.

Particle level simulations (PLS) have been used in the past for microstructure predictions in the plastics industry to learn about the real mechanisms present in processing discontinuous fibers [20,21,22]. PLS provide a fundamental modeling approach, as they solve the forces and moments acting on individual fibers during processing by accounting for the complex interactions of fibers and matrix. Compared to continuum models, the state and motion of fibers are not described as averaged and homogenized properties, but rather by solving the governing equations of each fiber explicitly. Hydrodynamic forces, fiber flexibility, excluded volume forces, fiber–fiber, and fiber–wall contacts are taken into account, leading to more accurate simulation results [7,23]. Additionally, PLS can be used to determine the fitting parameters of continuum models numerically. This has advantages as all parameters can be accurately controlled, detailed information is always available, and simulations are relatively inexpensive to perform.

PLS of semi-concentrated to concentrated fiber suspensions must address two central problems: First, the formulation of the equation of motion for individual fibers. Second, the particle–particle interactions, which can affect the flow behavior. To do this, various approaches have been taken in the past. Yamane et al. [24] represented fibers a rigid rods and calculated motion based on Jeffery’s theory. They considered only short range interactions and modeled them with lubrication theory. Fan et al. [20] extended Yamane’s work with rigid fibers. They accounted for long range interactions by employing a slender body approximation. Fan and co-workers later used a chain of beads joined with connectors to model flexible fibers, improving their viscosity predictions [25]. Schmid et al. [26] used PLS to study fiber flocculation. They modeled fibers as chains of rigid rods interconnected with hinges. Rods could rotate and twist about the hinges, replicating fiber bending and twisting deformations. They modeled interactions by considering repulsive forces acting normal to fiber surfaces to represent the fiber’s excluded volume, and frictional fiber forces which prevented fibers from sliding over one another.

The objective of this work is to generate reliable FO evolution data in a well-defined simple shear flow to aid in the validation and development of our PLS. Simple shear was chosen as it is one of the fundamental flow conditions present in most polymer processes, such as injection molding (IM). This allows us to directly correlate the rate of deformation with the filler’s behavior. IM and compression molded (CM) glass fiber-reinforced polypropylene samples were sheared in a Sliding Plate Rheometer (SPR) following Cieslinski et al. [27]. As has been shown by the same author, CM has no control over the planar orientation of the fibers, and is therefore not a suitable sample preparation method. This work will present a CM technique which ensures a controlled and repeatable initial FO for shear experiments. Results from both simulation and experiment are presented in this work.

2. Direct Fiber Model

The direct fiber or mechanistic model, used in this work, was developed at the Polymer Engineering Center, University of Wisconsin-Madison, USA, and is based on Schmid et al. [26]. Every single fiber is discretized as a chain of rigid elements with circular cross sections of diameter D, interconnected with spherical joints. The joints elastically couple the segments of a fiber as shown in Figure 1.

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⁻³. In our work, the Re was calculated to be in a range of 10⁻⁹–10⁻⁴. 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³ 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_i^H$, the interaction force with an adjacent segment $j$ $F_{ij}^C$, and intra-fiber forces $X_i$ exerted by internal connection loads at the nodes. The translation equation of motion is

$$0=F_i^H+{{\displaystyle \sum }}_j{F}_{ij}^C+X_i-X_{i+1}$$

(1)

The rotational equation of motion is equally derived and includes elastic recovery terms $M^b$ and a hydrodynamic torque $T_i^H$

$$0=T_i^H-r_i\times X_{i+1}+{{\displaystyle \sum }}_j{r}_{ij}\times F_{ij}^C+M_i^b-M_{i+1}^b$$

(2)

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

$$0=v_i+{\omega }_i\times \left(x_{i+1}-x_i\right)-v_{i+1}$$

(3)

Hydrodynamic effects are approximated by modeling each fiber segment as a chain of beads. The hydrodynamic force $F_i^H$ is the sum of forces experienced by the beads $F_k^H$ and is described as

$$F_i^H={\displaystyle \sum }_{k=1}^m{F}_k^H={\displaystyle \sum }_{k=1}^{m}6\pi \mu a\left(U_k^{\infty }-u_k\right)$$

(4)

$F_k^H$ is given by Stokes law, where $k$ describes the number of beads, $a$ the bead radius, $\mu$ the polymer matrix viscosity, $U_k^{\infty }$ the surrounding fluid velocity, and $u_k$ the velocity of the bead $k$ ($u_k=u_i+{\omega }_i\times r_k$).

$F_i^H$ is then be written as

$$F_i^H={\displaystyle \sum }_{k=1}^{m}6\pi \mu a{U}_k^{\infty }-m{u}_i-{\omega }_i\times {\displaystyle \sum }_{k=1}^m{r}_k$$

(5)

Due to the fluid’s vorticity each bead is also subjected to a hydrodynamic torque $T_i^H$

$$T_i^H={\displaystyle \sum }_{k=1}^m{T}_k^H$$

(6)

where the hydrodynamic contribution $T_k^H$ of bead $k$ is

$$T_k^H=8\pi \mu a^3\left({\Omega }_k^{\infty }-{\omega }_k\right)$$

(7)

${\Omega }_k^{\infty }$ is the vorticity of the surrounding fluid and ${\omega }_k$ is the angular velocity of the bead $k$. Substituting the expression of $T_k^H$ and the expression of $u_k$ into the fluid’s vorticity we can write

$$T_i^H={\displaystyle \sum }_{k=1}^m\left(8\pi \mu a^3{\Omega }_k^{\infty }\right)-8\pi \mu a^{3}m{\omega }_i+{\displaystyle \sum }_{k=1}^{m}6\pi \mu a\left(r_k\times U_k^{\infty }\right)+u_i\times {\displaystyle \sum }_{k=1}^m(6\pi \mu a{r}_k)-{\displaystyle \sum }_{k=1}^m(r_k\times ({\omega }_i\times r_k))$$

(8)

The fiber–fiber interaction force $F_{ij}^C$ is decomposed in a normal force $F_{ij}^N$ and a tangential force $F_{ij}^T$

$$F_{ij}^C=F_{ij}^N+F_{ij}^T \mathrm{for} d_{ij} \text{<} d_{threshold}$$

(9)

$F_{ij}^C$ 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_{ij}^N$ is an excluded volume force which is implemented as a discrete penalty method [7,32,33]. $F_{ij}^N$ usually has a function of the following type [26,34,35],

$$F_{ij}^N=A exp\left[-B\left(\frac{2d_{ij}}{D}-2\right)\right]n_{ij}$$

(10)

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_{ij}^T$ between segment $i$ and $j$ is computed as

$$F_{ij}^T={\mu }_f \left|F_{ij}^N\right| \frac{\mathsf{\Delta }{u}_{ij}}{\left|\mathsf{\Delta }{u}_{ij}\right|}$$

(11)

${\mu }_f$is the Coulomb coefficient between fibers, $\mathsf{\Delta }{u}_{ij}$ is the vector of the relative velocity between elements.

Bending of a fiber was approximated by using elastic beam theory

$$M_i^b=\frac{\left(\pi -\alpha \right)IE}{l}\left(\frac{r_i\times r_{i-1}}{\left|r_i\times r_{i-1}\right|}\right)$$

(12)

with $M_i^b$ as the bending moment, $\alpha$ 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].

3. Test Setup and Experimental Validation

3.1. 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³) are chemically coupled to the PP matrix (ρ = 0.91 g/cm³). 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. Fiber properties. Average FLs, a₁₁, a₂₂, and a₃₃ as the orientation tensors. a₁₁ representing extrusion direction in compression molded samples and the melt flow direction in injection molded samples.

Material Property	CM Plates	IM Plates
LN [mm]	0.83	1.28
LW [mm]	1.53	2.92
a11 [-]	0.86	0.60
a22 [-]	0.11	0.37
a33 [-]	0.03	0.03

.

The second sample preparation procedure involved IM a simple plate geometry (102 × 305 × 2.85 mm³). 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.

3.2. 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² 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 a₁₁ and the a₂₂ 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⁻¹. 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.

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

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

4. Direct Fiber Simulation

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.

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

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⁻¹ 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. Cluster properties for injection- and compression molded samples.

Cluster Property	IM Plates	CM Plates
vol %	4	8.5
Longest fiber [mm]	5	4
LN [mm]	0.71	0.83
LW [mm]	1.32	1.53
a11 [-]	0.6	0.36
a22 [-]	0.39	0.62
a33 [-]	0.002	0.02

.

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. Shear cell properties for injection- and compression molded samples.

Parameter	Value
E [GPa]	73
Fiber diameter [µm]	19
η [Pas]	110
$\dot{\gamma }$ [s−1]	1
Time step [s]	5 × 10−5
Integrations [-]	1,200,000

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

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 (a₂₂), while fibers in the shell are oriented along the flow direction (a₁₁) due to the fountain flow effect. The orientation in thickness direction a₃₃ 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 a₁₁ 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₁₁ 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₁₁ before fibers start aligning in shearing direction. This drop is less pronounced in the PLS, as fibers start off with a higher a₁₁ alignment. An initial decrease in a₁₁ 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₁₃, thus a decrease in a₁₁. This increase in a₁₃ could be seen in the experimental data. Depending on the initial value of a₁₃, the alignment in flow direction is different. If a₁₃ is positive, the fiber will continue to align in flow direction. However, an initial negative a₁₃ orientation would require the fiber to flip before becoming aligned in flow direction. It is expected that a fiber with a negative a₁₃ would take longer to align in flow direction [27]. A slight initial decrease in a₁₁ 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 a₁₁ 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.

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 a₁₁ could be observed (Figure 9). CM samples with a controlled FL, FC, homogenous structure, and with a low a₁₁ 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 a₁₁ 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₁₁ 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 $\kappa$ 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 $\left(1-\varphi \right)$) 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₁₃ and a₁₂ components of orientation are shown in Figure 12. In both prediction and experiment a₁₃ and a₁₂ show close agreement with the final steady state value. The measured a₁₃ component shows small oscillations. This behavior indicates some fibers orbit in the plane of shear. The predicted a₁₃ component also shows these orbits, however, with a shorter period. The initial measured a₁₃ component has a negative value of −6 × 10⁻⁴, while the computational cluster started with a very small, yet positive, value of 6 × 10⁻⁵. Literature suggests that the initial value of a₁₃ 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₁₁ component oriented directly to steady state in the simulation, while the experimental a₁₁ did show a slight undershoot (Figure 11). It can be seen that a₁₂, for the simulation, increases continuously with shear strain. As fibers oriented in the 1–2 plane orient themselves in 1, the value of a₁₂ 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 a₁₁ 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 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 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₁₁ 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 $\varphi r_p$.