**1. Introduction**

Molding of discontinuous reinforced polymer fiber suspensions, e.g., glass fibers in polymer melt, leads to fiber reorientation. Understanding and predicting the behavior of such fiber suspensions is crucial for achieving high quality products in common composite manufacturing processes such as injection molding and compression molding. Due to the high cost of production facilities and molds, it is desirable to simulate the flow in a computationally efficient and reliable way before running experiments. The fiber reorientation is also of high interest for subsequent anisotropic structural simulations in the framework of a continuous CAE chain [1,2].

Jeffery [3] was the first to analytically derive the motion of a single rigid, ellipsoidal shaped body in a viscous Newtonian flow without buoyancy or inertia. Jeffery's work was extended by Folgar and Tucker [4] to account for fiber interactions by introducing a fiber interaction coefficient that models isotropic diffusivity of fiber orientation. Other phenomenological parameters were introduced to capture experimentally observed orientation delays in the SRF and RSC model [5] and to account for anisotropic diffusion in the Anisotropic Rotary Diffusion (ARD) model [6] or the improved Anisotropic Rotary Diffusion (iARD) model [7]. These phenomenological parameters account for interactions of multiple fibers in non-dilute suspensions and are fitted to experimental observations, but do not describe a two-phase suspension. Instead, they model the fiber orientation state with second and

fourth order moments of the fiber orientation distribution function <sup>Ψ</sup>(**p**) introduced by Advani and Tucker [8] as fiber orientation tensors

$$\mathbf{A} = \int\_{S} \mathbf{p} \otimes \mathbf{p} \, \mathbb{1} \, \mathbf{p} \, (\mathbf{p}) \, \mathrm{d}p \tag{1}$$

and

$$\mathbb{A} = \int\_{S} \mathbf{p} \otimes \mathbf{p} \otimes \mathbf{p} \otimes \mathbf{p} \,\,\Psi(\mathbf{p}) \,\mathrm{d}p.\tag{2}$$

Here, **p** describes a fiber direction and d*p* is the surface element on a unit sphere *S* := {**p** ∈ R<sup>3</sup> : **p** = <sup>1</sup>}. Typically, a closure approach for the fourth order fiber orientation tensor is employed to describe the evolution of the second order fiber orientation tensor. Whenever this work utilizes a closure approximation, the invariant-based optimal fitting (IBOF) closure is chosen [9]. Determining the parameters in these macroscopic models from experiments can be time consuming. Thus, a direct fiber simulation may be utilized to identify these parameters.

#### *1.1. Point-Wise Interaction Methods*

A common approach for the simulation of fiber suspensions is the treatment of fibers as slender bodies that interact with the fluid at discrete points. Exact solutions from Stokesian dynamics [10] or slender-body theory [11] are utilized to describe long-range hydrodynamic interaction between fluid and particles for creeping flows. Several authors proposed models for single flexible fibers suspended in a fluid. Hinch [12] started by modeling inextensible threads. Yamamoto et al. [13,14] developed a fiber model consisting of individual beads that experience Stokesian drag forces from the fluid. These bead chain models are computationally expensive and authors have combined multiple beads to rods leading to more efficient rod chains [15–17]. Alternatively, spheres [18,19] and spheroids [20] connected with hinges were suggested. Lindström and Uesaka [21] use a discrete field of point forces to ensure that the fluid is experiencing forces opposed to the fiber (two-way coupling). Several authors investigated rheological properties of multiple rigid fibers suspended in a fluid based on these models for single fibers. Yamane et al. [22] described the motion of multiple rigid rods with hydrodynamic drag force and torque on the individual rods. They added a lubrication force for rods that are in close proximity to each other in order to capture short range hydrodynamic effects. However, they did not account for long range hydrodynamic interaction between particles, which was addressed later by Fan et al. [23] using slender-body theory. Sundararajakumar and Koch [24] showed that lubrication forces alone do not prevent penetration of fibers and included contact forces. Suspensions of multiple flexible fibers were investigated with rod-chain models [25] and simplified bead chain models [26]. In addition to these rheological investigations, few direct simulations have been applied on component scale with flexible fibers or bundles to investigate effects such as fiber-matrix separation [27–30].

#### *1.2. Resolved Methods*

In contrast to models based on discrete point-wise interaction forces, the suspension flow might be also described directly by a two-phase flow, in which one phase represents fibers and the other one the suspending fluid. Solving this fully coupled two-phase system with mesh-based approaches such as Finite Elements or Finite Volumes raises the problem of large mesh deformations. A simulation using an immersed boundary method resolves the remeshing problem but requires a mesh resolution significantly smaller than the fiber diameter to track the interface and is therefore computationally expensive. An alternative approach for the simulation of two-phase fiber suspensions is the use of particle methods. Meshless particle methods can offer a natural way to represent fluid-structure interaction at the fiber surface and have been studied less than point-wise interaction methods. Bian and Ellero propose a splitting scheme for separate integration of long-range hydrodynamic interactions and short-range lubrication [31], which was applied in an SPH simulation of a concentrated spherical

particle suspension [32]. The fine resolution of suspended particles leads to a high accuracy, but comes at high computational costs. Yashiro et al. [33,34] connected particles to rigid bodies for modeling the injection molding of dilute short-fiber-reinforced composites using a moving particle semi-implicit method. Recently, a very similar approach using SPH was reported by He et al. [35,36]. However, both showed only short time periods when comparing the simulation to Jeffery's equation and did not investigate the interaction between fibers. Fiber suspensions were also investigated in the framework of dissipative particle dynamics with a focus on Boger fluids [37]. This investigation used a Langrange multiplier to constrain fiber extension and multiple parameters had to be calibrated to analytical solutions in order to achieve correct hydrodynamic forces on the fiber. Yang et al. [38] employed the SPH method in order to simulate a single flexible fiber in a viscous fluid with a focus on determining its drag coefficient. However, they limited their work to a single fiber without interactions.

The present work extends Yang's model with a viscous surface traction term necessary to meet Jeffery's equation precisely and allows for the simulation of non-dilute fiber suspensions through the introduction of fiber contact forces. First, rotation and bending of a single fiber is compared to Jeffery's equation to validate the approach. Then, parameters of the RSC fiber orientation model are determined using SPH simulations.
