**2. Theory**

The orientation of a single fiber can be characterized by a unit vector p (see Figure 1 [4]).

$$p\_1 = \sin \Theta \cos \phi \tag{1}$$

$$p\_2 = \sin\Theta \sin\phi \tag{2}$$

$$p\_3 = \cos \Theta.\tag{3}$$

**Figure 1.** Definition of a single unit fiber in space.

Various models have been developed to model fiber orientation. Work was performed at Virginia Tech [5] using a sliding plate rheometer to choose a fiber orientation model—the pARD-RSC (principal anisotropic rotary di ffusion-reduced strain closure) [6] model. The chosen model was implemented in AMI ® using the Solver API solution within AMI ® [5]. The model is given by

.

$$\begin{split} \frac{\partial \mathbf{A}}{\partial \mathbf{I}} = \dot{A} = \mathbf{W} \cdot \mathbf{A} \quad -\mathbf{A} \cdot \mathbf{W} + \zeta [\mathbf{D} \cdot \mathbf{A} + \mathbf{A} \cdot \mathbf{D} - 2(\mathbf{A} + (1 - \kappa)(\mathbf{L} - \mathbf{M} : \mathbf{A})) : \mathbf{D}] \\ \qquad + \dot{\chi} [2[\mathbf{C} - (1 - \kappa)\mathbf{M} : \mathbf{C}] - 2\kappa(tr \mathbf{C})\mathbf{A} - \mathbf{5}(\mathbf{C} \cdot \mathbf{A} + \mathbf{A} \cdot \mathbf{C}) \\ \qquad + 10[\mathbf{A} - (1 - \kappa)(\mathbf{L} - \mathbf{M} : \mathbf{A}) : \mathbf{C}] \end{split} \tag{4}$$

$$\begin{aligned} \mathcal{L} &= \sum\_{i=1}^{3} \lambda\_i \mathbf{e}\_i \mathbf{e}\_i \mathbf{e}\_i \\ \mathcal{M} &= \sum\_{i=1}^{3} \mathbf{e}\_i \mathbf{e}\_i \mathbf{e}\_i \mathbf{e}\_i \\ \mathcal{C} &= C\_1 \mathcal{R}\_A \left( \begin{array}{cccc} D\_1 & 0 & 0 \\ 0 & D\_2 & 0 \\ 0 & 0 & D\_3 \end{array} \right); \text{ with } A = \mathcal{R}\_A \left( \begin{array}{cccc} \lambda\_1 & 0 & 0 \\ 0 & \lambda\_2 & 0 \\ 0 & 0 & \lambda\_3 \end{array} \right) \end{aligned} \tag{5}$$

$$D\_1 = 1, \; D\_2 = \mathbf{c}, \; D\_3 = 1 - \mathbf{c}$$

where *A* is the second order orientation tensor, and A is the fourth order orientation tensor; *RA* is the Eigen matrix, *D* is the strain rate, and W is the vorticity tensor. *CI* is the interaction coefficient introduced in [7].

The parameters *CI*, *D*2, and κ were determined from a numerical shear cell experiment developed at the Polymer Engineering Center (PEC) of the University of Wisconsin—Madison [8]. For the generation of unit cells with discrete fibers, sequential addition and the migration method were used [9]. A statistical approach to determine the uncertainties of this procedure was established by varying the initial fiber cell parameters.

#### **3. Generation of the Unit Cells**

Two different setups of cells were generated. First, the influence on cell generation was investigated. The generation of the cells containing 30 wt%. in a 512 μm × 512 μm × 512 μm cubic cell resulted in 1118 discrete fibers with a diameter of 10 μm and a fiber length of 275 μm. The algorithm used assures that the fibers do not bend or overlap with each other [9]. The initial orientations, given by the fiber orientation tensor elements *A*11, *A*22, and *A*33, were 0.15, 0.05, and 0.8 and correspond to the initial alignment of pellets in the sliding plate experiments. All the diagonal elements, such as *A*12, were set to 0. Under these settings, 10 initial cells were created by the fiber generator. Notably, the second order orientation tensor does not impose a unique orientation state. Consequently, the orientation between the 10 cells was not identical.

Second, the influence of fluctuations in the initial settings for orientation, fiber length, and fiber content were examined. A Design of Experiment (DoE) scheme was used to investigate three parameters (*A*33, fiber length, and fiber concentration) in three levels. To minimize the calculation times, a reduced DoE scheme (D-optimal using the software Cornerstone®) was applied, resulting in 15 settings to consider linear and quadratic combinations and interactions. The upper and lower bounds were chosen as the upper and lower bounds of the material properties of PBT GF 30 during injection molding, since this study aims to determine the influence of varying the parameters during the final fiber orientation of the part.

The settings applied for the D-optimal scheme with three parameters and three steps are given in Table 1.



Figure 2 shows the generated fibers for V15 from Table 1 prior to the applied shear.

**Figure 2.** Initial fiber alignment for V15.

The mechanistic modelling approach of PEC was applied to study the orientation of the fibers in a shear flow field. The numerical shear cell applies shear in the 1, 2 plane. On the top and bottom of the shear cell, walls that prevent fibers from penetrating in two directions are modelled. For the 1 and 3-directions, periodic boundary conditions were applied to the cell. A shear rate of 1 s<sup>−</sup><sup>1</sup> was applied for 100 s and 200 s. The input parameters for the calculations include a constant viscosity (taken for a fixed temperature with the shear rate from the data sheets for an unreinforced PBT), the fiber geometry, the fiber volume content, and the cell dimensions, as well as the number of elements. The fibers are modelled as stiff rigid bodies with a single element.

In each time step, the velocity, angular velocity, and position of each fiber was calculated with the force and momentum balance. The forces/moments considered in the model are hydrodynamic forces exerted from the polymer on the fibers and interaction forces between fibers. Back coupling between fiber motion and fluid motion was not considered.

Between two approaching fibers, it is assumed that the interaction is divided into three regimes: lubrication forces at large distances, transition forces at the distance of the surface roughness of the fibers, and mechanical forces at contact [10]. The lubrication forces are dependent on the distance and approaching angle. The lubrication force model is based on the work in [10]. Since analytical solutions exist only in the parallel case [11] and in the case of infinitely long fibers [12], optimization of the lubrication force was conducted. As an experimental approach was not available, a numerical one was applied. A fully coupled simulation in COMSOL® was used to calculate the lubrication force between two fibers at varying angles, relative velocities, fiber lengths, and matrix viscosities. Based on these data, an analytic optimization of the lubrication force was performed and implemented in the mechanistic model.

#### **4. Results and Discussion**

After applying a constant shear rate of 1 1/s for times of 100 s or 200 s to the initially aligned shear cell, the evolution parameters according to Equation (4) were determined based on a generic fitting algorithm in MATLAB®. Figure 3 shows a graph of the fiber orientation evolution over time for V15. This graph shows the evolution of fiber orientations for *A*11, *A*22, and *A*33 over strain under a constant shear rate. The parameters for *CI*, *<sup>D</sup>*2, and κ are derived from the fitted lines according to Equation (4).

**Figure 3.** Fiber orientation evolution for a shear rate of 1 s<sup>−</sup><sup>1</sup> and parameters set as V15.

From the variation of cell generation with identical input parameters, it can be found that *CI* varies significantly with the generation of cells and model parameter fitting. The mean values determined over all ten generations can be found together with the standard deviation in Figure 4.

**Figure 4.** Average value for ten cells with identical inputs.

To determine whether such deviations play an important role, DoE calculations were performed. Table 2 shows the resulting parameters for each run.



A statistical evaluation of the three factors (F1 = length, F2 = content, F3 = orientation) and the three response data (R1 = *CI*, R2 = *D*2, R3 = κ) was performed and is plotted in Figure 5a,b. For the chosen parameter range, no clear correlation could be found. The response (R1 to R3) was independent from the input parameters (factors F1 to F3).

**Figure 5.** (**a**) Predicted Response Graph for R1 to R3 vs. F1 to F3. (**b**) Influence of factors on the response *CI* (top), *D*2 (middle), and κ (bottom).

In the Pareto graphs, for *CI* and *D*2, the fiber content has the largest positive effect, and κ is the initial orientation. Nevertheless, the effect, in general, is rather small.

#### **5. Numerical Simulation of the Fiber Orientation in Moldflow**

The influence of these varying parameters on the calculated fiber orientation in AMI® 2019 was investigated with seven parameter sets from V3, V6, V8, V9, V10, V13, and V14. These seven settings were chosen because they represent the extreme values from all valid parameter sets (min/max).

Simulations were done using the Solver API functionality in AMI® 2019 in a pre-release version. The geometry used for these simulations was a center-gated box with a non-homogeneous wall thickness distribution (see Figure 6).

**Figure 6.** Geometry of the simulated center gated box.

For the simulation, material data from the AMI® 2019 data base for Ultradur B4300 G6 were chosen. The injection parameters used are given in Table 3.



The cooling system of the mold was also taken into account (see the blue cylinders in Figure 7).

**Figure 7.** Model in AMI® 2019 with the cooling system (blue) and melt entrance (yellow cone).

To make a quantitative comparison between the seven chosen model parameter sets, each was simulated with identical processing parameters given in Table 3. Two element stacks of hexaeder elements over thickness were generated. The calculated fiber orientations were then mapped on these stacks. The positions of each stack are given in Figure 8 (identified by red circles).

**Figure 8.** Positions of the element stacks to evaluate fiber orientation.

For each calculation, the deviation to the first (V3) model parameter set was determined according to the root mean square deviation (RMSD) method.

$$\text{RMSD} = \sqrt{\frac{\sum\_{n=1}^{N} \left(\chi\_{1,n} - \chi\_{2,n}\right)^2}{N}} \tag{6}$$

Table 4 gives a summary of the RMSD determined for each position.



With a maximum orientation of 0.8 in the flow direction, for the chosen input parameters and the determined model parameters of the DoE, the change in fiber orientation is 10–20%.
