**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–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–14], length distribution [15,16], and fiber content prediction [17–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–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 a ffect 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 Je ffery'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.

**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<sup>3</sup> 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 segmen<sup>t</sup> *i* includes the drag force *FHi* , the interaction force with an adjacent segmen<sup>t</sup> *j FCij*, and intra-fiber forces *Xi* exerted by internal connection loads at the nodes. The translation equation of motion is

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

The rotational equation of motion is equally derived and includes elastic recovery terms *M<sup>b</sup>* and a hydrodynamic torque *THi*

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

where *rij* 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 + a\_i \times (\mathbf{x}\_{i+1} - \mathbf{x}\_i) - v\_{i+1} \tag{3}$$

Hydrodynamic effects are approximated by modeling each fiber segmen<sup>t</sup> as a chain of beads. The hydrodynamic force *FHi*is the sum of forces experienced by the beads *FHk*and is described as

$$F\_i^H = \sum\_{k=1}^m F\_k^H = \sum\_{k=1}^m 6\pi\mu a \left(\mathcal{U}\_k^{\infty} - u\_k\right) \tag{4}$$

*FHk* is given by Stokes law, where *k* describes the number of beads, *a* the bead radius, μ the polymer matrix viscosity, *<sup>U</sup>*<sup>∞</sup>*k*the surrounding fluid velocity, and *uk* the velocity of the bead *k* (*uk* = *ui* + ω*i* × *rk*). *FHi*is then be written as

$$F\_i^H = \sum\_{k=1}^m \epsilon \pi \mu a L\_k^{\infty} - m\mu\_i - \omega\_i \times \sum\_{k=1}^m r\_k \tag{5}$$

**Figure 2.** Fiber represented as a chain of segments. Positions are stored at the segmen<sup>t</sup> nodes. A neighboring fiber is depicted. The contact force *FCij* - 0 when the distance between two segments is *dij* < *D* [32].

Due to the fluid's vorticity each bead is also subjected to a hydrodynamic torque *THi*

$$T\_i^H = \sum\_{k=1}^{m} T\_k^H \tag{6}$$

where the hydrodynamic contribution *THk*of bead *k* is

$$T\_k^H = 8\pi\mu a^\beta \left(\Omega\_k^{\alpha \circ} - a\_k\right) \tag{7}$$

<sup>Ω</sup><sup>∞</sup>*k* is the vorticity of the surrounding fluid and ω*k* is the angular velocity of the bead *k*. Substituting the expression of *THk*and the expression of *uk* into the fluid's vorticity we can write

$$T\_i^H = \sum\_{k=1}^m \left( 8\pi\mu a^3 \Omega\_k^{\prime \alpha} \right) - 8\pi\mu a^3 m\omega\_i + \sum\_{k=1}^m 6\pi\mu a \left( r\_k \times \mathcal{U}\_k^{\prime \alpha} \right) + \mu\_i \times \sum\_{k=1}^m \left( 6\pi\mu a r\_k \right) - \sum\_{k=1}^m \left( r\_k \times \left( \omega\_i \times r\_k \right) \right) \tag{8}$$

The fiber–fiber interaction force *FCij*is decomposed in a normal force *FNij*and a tangential force *FTij*

$$F\_{ij}^C = F\_{ij}^N + F\_{ij}^T \text{ for } d\_{ij} < d\_{threshold} \tag{9}$$

*FCij* describes the friction between the segments and starts acting when the distance between adjacent fibers, *dij*, is below a defined threshold *dthreshold*. *FNij* is an excluded volume force which is implemented as a discrete penalty method [7,32,33]. *FNij* 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} \tag{10}$$

where *dij* is the shortest distance between rod *i* and rod *j*, *D* the fiber diameter, and *nij* 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 *FTij*between segmen<sup>t</sup> *i* and *j* is computed as

$$F\_{ij}^T = \mu\_f \left| F\_{ij}^N \right| \frac{\Delta u\_{ij}}{\left| \Delta u\_{ij} \right|} \tag{11}$$

μ*f* is the Coulomb coefficient between fibers, <sup>Δ</sup>*uij* is the vector of the relative velocity between elements.

Bending of a fiber was approximated by using elastic beam theory

$$M\_i^b = \frac{(\pi - \alpha)IE}{l} \left(\frac{r\_i \times r\_{i-1}}{|r\_i \times r\_{i-1}|}\right) \tag{12}$$

with *Mbi*as the bending moment, α as the angle between segments, and *l* as the segmen<sup>t</sup> 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 SABICTM. 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/cm3) are chemically coupled to the PP matrix ( = 0.91 g/cm3). 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 (SABICTM 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 (LN = 6.7 mm, LW = 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 LN of 0.83 mm and LW 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.



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

**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, e ffectively 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<sup>−</sup>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.

#### *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<sup>3</sup> 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.

**Figure 4.** Fully characterized microstructure of an injection molded plaque.

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